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ABSTRACT 


The potential of a relatively new method to remotely measure near-surface cur¬ 
rents and dominant wave spectra using Interferometric Synthetic Aperture Radar 
(INSAR) is demonstrated. INSAR consists of a single conventional SAR augmented 
by an additional receiving antenna. The phase difference between the correspond¬ 
ing SAR image scenes observed by the antennas provides an interferogram that is 
directly proportional to the ocean surface velocity field. This direct motion de¬ 
tection by INSAR suggests a significant advance compared with conventional SAR 
where the response to the moving ocean surface is indirectly related to the complex 
modulation of the surface reflectivity by longer waves and currents. 

An experiment using an airborne implementation of INSAR to measure ocean 
surface currents and wave fields, compared with simultaneous ground truth measure¬ 
ments using Lagrangian drifters and wave array data was conducted in Monterey 
Bay. INSAR measured mean current magnitude estimates agree to within 10 per¬ 
cent compared with conventional measurements. The INSAR image wavenumber 
spectrum is consistent with the in situ directional spectrum and with predicted 
numerical refraction model outputs. The wavelength of the observed swells are in 
better agreement (correlation better than 0.9) than wave direction. 

An attempt to estimate the scene coherence time for L-band SAR was made by 
taking advantage of the almost simultaneously acquired SAR and INSAR images. 
The obtained mean scene coherence time 0 (100 msec) is consistent with sparse 
observed estimates in the literature. This limited radar temporal coherence, caused 
by the velocity spread of short waves, degrades the azimuthal resolution of SAR 
and INSAR and depends on the sea state and radar wavelength. The experimental 
results show that the finite scene coherence time has a dominant role on the distor¬ 
tion of the INSAR (and SAR) image spectrum relative to the ocean wave spectrum. 
The present study introduces limitations c! interferometric SAR configuration in 
imaging nonstationary scenes like the ocean surface. 
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I. INTRODUCTION 


A. INTERFEROMETRIC RADAR REMOTE SENSING 
OF OCEAN SURFACE VELOCITY FIELDS 

Historically, oceanography has depended on sparse observations made from 
ships. Now, orbiting platforms offer global coverage of the world’s oceans using 
various types of sensors in a broad spectrum (visible, IR, microwave). “How relevant 
are these surface observations in understanding today’s oceanographic problems? 
How well matched to the oceanic time and space scales are the satellite samples? 
What oceanic processes can be observed, and in what way? And, finally, what 
instruments and systems exist to provide oceanic data from space?” (Stewart, 1985). 

Some of the answers to these questions are addressed in this study by ap¬ 
plying a new method of remotely measuring ocean surface velocity fields, using a 
conventional Synthetic Aperture Radar (SAR) augmented with two spatially sep¬ 
arated antennas (Goldstein and Zebker, 1987). By incorporating the physics of 
interferometry, this unique remote sensing method offers an opportunity for esti¬ 
mating synoptic spatial structures of ocean near surface currents and wavefields with 
high grid nominal resolution (12 m in this work) and large coverage, never before 
obtained. 

The objectives of this research are to develop a theoretical foundation for 
Interferometric SAR (INSAR) in imaging the free surface, distributed moving ocean 
surface, to verify its feasibility with simultaneous ground truth experiments, and to 
develop interferometric SAR image processing techniques to estimate spatial synop¬ 
tic currents and wave spectra. An experiment was conducted to ground truth the 
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velocity field imagery of INSAR using surface Lagrangian drifters, wind measure¬ 
ments, and a shallow water wave array. 

B. STATEMENT OF OBJECTIVES AND PROBLEMS 

Three specific objectives are outlined for this dissertation: 

1. To develop a theoretical basis for the capability of SAR in an interfero¬ 
metric mode to measure the ocean surface velocity field. The basic limitations in 
correlating the backscatter radar echo with ocean phenomena is that the ocean sur¬ 
face combines random fluctuating, non-stationary and non-homogeneous scattering 
elements. In this study, a preliminary theoretical model is developed which repre¬ 
sents the major physical parameters involved in remotely measuring near surface 
currents and wavefields using synthetic aperture radar in an interferometric mode. 

The decorrelation of the ocean surface scene (finite scene coherence time) 
as a function of INSAR parameters is quantified. This allows to characterize the 
sensitivity and capability of INSAR in terms of imaging near surface ocean currents 
and wave fields. Noncoherent motion effects on airborne SAR (and hence, INSAR) 
imagery, which is an area of controversy, is also addressed in the present study. 

2. To develop analysis techniques to obtain quantitative spatial information 
of near ocean surface currents and wave fields from INSAR imagery. Unlike SAR, 
which has serious shortcomings in obtaining directional wave spectra, INSAR has 
an inherent improved potential of measuring 2D wavenumber spectrum of the ocean 
surface dominant wave field. This improvement is related to the direct motion de¬ 
tection mechanism of INSAR when sufficient Bragg scatterers are present. In this 
study, appropriate theoretical tools are developed to separate wave field velocities 
from ocean currents. This separation allows accurate and unique information con¬ 
cerning spatial structures of dominant ocean wavefields, currents, and other features. 
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3. To ground truth the currents and wave fields measured by INSAR. A si¬ 
multaneous in situ and radar remote sensing experiment was conducted to determine 
the feasibility of INSAR to image near surface ocean dynamics compared with more 
conventional measurements. The analyzed data also provide a unique study of wave 
refraction and shoaling and currents in the complex nearshore region of Marina in 
Monterey Bay. 

C. SUMMARY 

The present study demonstrates the feasibility of this unique technique to 
provide reliable quantitative spatial information about dynamic processes at the 
ocean surface. Realization of the potential of INSAR with spacecraft might make 
it possible to routinely monitor global ocean wave spectra and surface currents. A 
review of electromagnetic interaction with oceanic backscattering waves is presented 
in Chapter II along with backscattering models which supports the interpretation 
of the results. The principal mechanisms of SAR and INSAR and relevant theory 
in imaging ocean surface dynamics are presented in Chapter III. In Chapter IV, a 
preliminary theoretical analysis for INSAR in imaging ocean surface velocity fields is 
presented for the first time. Based on the SAR two-scale model, a theoretical model 
for estimating the temporal decorrelation of the ocean surface using simultaneous 
SAR and INSAR observations is also developed. Chapter V describes oceanographic 
relevant theory which is essential for the interpretation of the results. A numerical 
refraction model (Dobson, 1967) is used to predict wave length and direction in 
shallow water. These outputs are compared with INSAR spectral analysis estimates 
(Chapter VII). Two-dimensional wavenumber spectra for in situ and interferometric 
SAR are successively described. A description of the experiment, environmental 
conditions, and data processing follows in Chapter VI. Data analysis methods and 
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experimental results are presented in Chapter VII. Results, analysis, and discussions 
are given in Chapter VIII. These findings are summerized and the major conclusions 
and suggestions for further study are stated in Chapter IX. 

J 
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II. LITERATURE REVIEW 


A. INTRODUCTION 

Due to rapid advances in microwave remote sensing of ocean surface during 
the last decade, much progress has been made in understanding the physics of 
scattering of microwave radiation from the ocean surface (Ilasselmann et al., 1985; 
Thompson, 1989). The result of this progress yields a capability to observe in more 
detail near-surface ocean features such as: wave height, wind speed, spatial wave 
field data, tides, eddies and currents (Valenzuela, 1989). 

The interpretation of microwave remote sensing data is done through the 
physics of microwave scattering theories. Therefore, the discussion in this chapter 
begins with the mechanisms of electromagnetic scattering from ocean surface for 
real aperture radars, which is the basis of radar scattering theory. Next, the dis¬ 
cussion extends to coherent synthetic aperture imaging radars, which in turn adds 
the comple motion effects when sensing the ocean. The complexity is mainly a 
result of the random uncorrelated motion of individual scatterers, smaller than the 
spatial resolution cell of SAR and of the order of the radar wavelength. The inter¬ 
action between different scales of waves and currents smear and defocus the SAR 
image, resulting in difficulties in correlating the image with the sensed phenomena 
(see Chapter IV). In spite of the difficulties, considerable progress has been made in 
understanding SAR imagery of ocean surface features. (Raney, 1971; Elachi, 1977; 
Alpers and Rul'enach, 1979; Valenzuela, 1980; Plant and Keller, 1983; Tucker, 1984; 
Ilasselmann et al., 1985; Kasilingham and Shemdin, 1988). 
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B. ELECTROMAGNETIC FIELDS SCATTERED 
FROM THE OCEAN SURFACE 

The interaction of electromagnetic waves with ocean surface is a complex 
boundary value problem, for which exact solutions are not available to the present 
time (Barrick, 1968; Valenzuela, 1971). This is because the ocean surface is a com¬ 
plex distributed scatterer, varying in time and space, and also weakly nonlinear 
behaving as a random rough surface. In the linear approximation, the statistics of 
the ocean surface are homogeneous, stationary, and Gaussian. With these assump¬ 
tions, several fundamental electromagnetic backscatter models, in good agreement 
with observations, have been derived (Barrick and Peake, 1968; Barrick, 1968, 1972; 
Hasselmann, 1971; Kirshen, 1971; Daley, 1973; Valenzuela, 1974, 1978; among oth¬ 
ers). These models are based upon: 

• The high frequency methods, including geometric and physical optics which 
yield specular point models for smooth and undulating surface (Kerr, 1951). 

• The low frequency method, using the perturbation method for a slightly rough 
surface yields the Bragg scattering model (Rice, 1951). 

• Combination of high and low frequency methods renders the two-scale com¬ 
posite model for rough surfaces. 

The purpose of these fundamental models of microwave scattering from the 
ocean surface is to convert the received electric signal at the sensor into a mean¬ 
ingful geophysical signal. These models (specular, Bragg-Rice, and composite) are 
described extensively in a review paper by Valenzuela (1978), which remain the ba¬ 
sis for studies of ocean surface dynamics sensed by microwave sensors even in the 
present time (see Figure 2.1). Recently, the composite model has been extended to 
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incorporate a third intermediate scale utilized specifically for SAR imagery appli¬ 
cations. This theory is summarized in great detail in a consensus paper by Has- 
selmann, Raney, Plant, Alpers. Shuchmann, Lyzenga, Rufenach, and Tucker, 1985. 
The scattering of electromagnetic fields from random rough surface differs from de¬ 
terministic surfaces scattering by being statistical in nature. Mean values of the 
fields and power may be determined by ensemble averaging these quantities over 
the assumed or known statistics of the rough surface (Beckmann and Spizzichino, 
1963). 


Incident 

Wave 



Incident 

Wave 



Undulating surface 


Incident 

Wave 


Scattering 

Pattern 





Rough surface 


Figure 2.1: Illustration of EM scattering from ocean surface. 


Following Valenzuela (1978), the scattered power density is given by the mean 
magnitude of the Poynting vector, which is a cross product of the electric and 
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magnetic fields 


p. = \R'{ e. x h;) (2.i) 

where H* is the complex conjugate of the scattered magnetic field of the EM wave, 
and E, is the scattered electric field of the EM wave. The brackets denote ensemble 
average. 

A measure of scattering power from a rough surface can be related to the 
radar cross-section, which is proportional to the ratio of the scattered power from 
the surface to the incident power. For the far field: 

<r = Jim (4xfl*<lB.|*)/l«P) (2-2) 

where \Ei\ 2 and (|E,| 2 ) are the incident and mean scattered power respectively. R 
is the distance of the observation point from the scattering rough surface. 

In the following sections the specular, Bragg-Rice, and the composite scat¬ 
tering models are presented. 


to radar 



Figure 2.2: Specular point backscatter (from Stewart, 1985). 

1. Specular Point/Physical Optics Model 

In this model, the reflections are from surface slopes which are oriented 
normal to the line of sight of the radar (Hagfors, 1964). The principal assumptions 
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of this model are that the surface radius of curvature is almost everywhere much 
larger than the radar wavelength (physical optics or Kirchhoff method), and that 
multiscattering and/or shadowing are negligible (see Figure 2.2). The range of 
validity of physical optics, or the tangent plan approximation, was shown by Shmelev 
(1972) as: 

2 rkji cos 7 > 1 and ~rkn cos 3 7 » 1 (2.3) 

where r is the local radius of curvature of the surface, 7 and kn are the incidence 
angle and the wavenumber of the electromagnetic wave in free space respectively. 

The backscatter power from a single facet (local tangent plan) is a func¬ 
tion of the probability that the tangent plane/facet will be oriented perpendicularly 
to the incidence direction of the electromagnetic wave. Summing up the contribu¬ 
tions from all properly oriented facets within the resolution cell of the radar results 
in the total backscatter specular power. 

Applying the physical optics approximation, Barrick (1968) derived a 
specular model to predict the normalized radar cross section for backscatter from a 
finitely conductive rough surface: 


a„ = (ir/cos < 7)p(i)i,i),,)l/’(0)P 


(2.4) 


where p(r/ r ,»/y) is the joint-probability density of orthogonal slopes 7/ r and tj v of the 
sea surface 77 with respect to the incident EM wave and p(0) is the Fresnel reflection 
coefficient for normal incidence. The slopes are the horizontal derivatives of the 
local vertical displacements of water surface from mean level. 

A simplified form of Equation (2.4) is achieved by assuming that the 
ocean surface roughness has a Gaussian random distribution. Thus, for an upwind 
looking radar, Equation (2.4) becomes: 

a ° = 2S"S~ se ° 4 ^) exp l" tan2 ( 2 - 5 ) 
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where Si and Si are the mean square slopes of upwind and crosswind respectively. 
Verified by observations, Equation (2.5) seems to be a good approximation for near 
normal backscatter from an undulating surface. 

Radar observations for the ocean indicate that the specular backscatter 
is 2 to 3 times larger from a trough than from a crest (Yaplee, 1971; Valenzuela, 
1978). This could be due to “compressibility” of the short gravity-capillary waves 
at the crests of the long waves, resulting in a diffuse scatter which reduces the 
specular backscatter. The reflection properties of a surface are known to be reduced 
by increase in roughness. Away from normal incidence,' the backscatter is no longer 
specular, but diffuse, being produced by shorter waves which have steeper slopes 
(roughness). 

2. The Bragg-Rice Model 

For a slightly rough surface, to first order, the backscatter energy, in 
directions of at least 20 degrees from normal incidence, is proportional to the ripple 
•variance spectrum evaluated at the Bragg resonant wave number (Hasselmann et 
al., 1985) 

k = 2knsmy ( 2 . 6 ) 

where k is the ocean wavenumber, and kn sin 7 is the projection of the incident radar 
wave number onto the ocean surface. This diffraction grating (resonant) backscat- 
tering theory for e.m. waves, was first developed by Rice (1951) and demonstrated 
experimentally by Crombie (1954) with an IIF Doppler radar. The Doppler spec¬ 
trum showed two peaks shifted from the EM frequency by the positive and negative 
Bragg waves frequency, for waves propagating towards and away from the sensor. 

Since the resonant (Bragg-Rice) backscattering mechanism has been 
shown to be the dominant backscatter for the ocean surface for incidence angles 


10 





ranging from 20 to 70 degrees, typical for SAR remote sensing, this backscattering 
mechanism is reviewed in more detail, together with the composite surface model 
introduced in the next section. 

The main assumptions in deriving the first order Bragg-Rice model are: 

# Local matching of the fields at the surface with linearized boundary conditions. 
Hence, the surface vertical perturbations, h, of the Bragg ripples are small with 
respect to the radar wavelength (Rayleigh condition). 

# The number of resonant scatterers is large enough to allow a proper statistical 
description of the backscattered reflectivity. 

# The slopes of the Bragg waves are assumed small. Jordan and Lang (1972) 
argue that this condition represents the real ocean because water waves break 
long before they become large enough to violate the linear boundary condition. 


Based on these assumptions, the resonant scatter theory yields a nor¬ 
malized radar cross section (Wright, 1966; Barrick, 1972; Jordan and Lang, 1979): 

a 0 = 87r^|</ o (7)| 2 V>(2A:/{sin7,0) (2.7) 


where < 7 ,^( 7 ) are the first order scattering coefficients and ip(k x , k y ) is the two dimen¬ 
sional ocean wave number spectrum properly normalized. The indices, (i,j) denote 
transmitted and received polarizations respectively: 


9ij{ 7) = 


9 hh { 7) 
9vu{ 7) 


9 hv ( 7) 

gwi'y) 


( 2 . 8 ) 


Depolarization effects for the ocean surface are generally smaller (of higher order) 
than the copolarized contribution, therefore only scattering coefficients with un¬ 
changed polarization with respect to the transmitted polarization are considered 
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here: 


, ^ (e r -l)cos 2 7 

9hh{ 7 ) = 7---— 7775 

[cos 7 + (e r — 5111 * 7 ) 2 ] 

, . (s r - l)[e r (l + sin 2 7 ) - sin 2 7 ] cos 2 7 

9vv{ 7) - r . 2 m2 

[e r cos 7 + (e r — sm 7 ) 2 ] 

where e r is the relative complex dielectric constant. 


9vv{ 7 )= 


( 2 . 10 ) 


Away from normal incidence, Bragg resonant scattering yields an in¬ 
creasing polarization ratio with increasing incidence angle, with a larger vertical po¬ 
larization scattering coefficient. For a surface with perfect conductivity, e r = —zoo, 
which is a reasonable assumption for the ocean surface for lower radar frequencies 
(i.e., HF), Equations (2.9) and (2.10) simplify (Stewart, 1985): 


ShhW ~ cos 1 0 


( 2 . 11 ) 


9vv{0) ~ 1 + sin 2 0 (2.12) 

Therefore, for increasing incidence angles, the horizontal polarization Bragg scat¬ 
tering coefficient decreases in magnitude, and is more sensitive to incidence angles 
relative to vertical polarization scattering coefficient except for near grazing inci¬ 
dence where the trapped waves and shadowing become dominant. 

The backscatter resonant cross-section is directly proportional to the 
two-dimensional ocean wave number spectrum satisfying the Bragg resonant condi¬ 
tion (Equation 2.7), which in turn suggests the usefulness of microwave backscatter 
for probing ocean waves. If the Bragg water waves are not the dominant waves of the 
ocean, higher order Doppler shifted backscatter will occur whenever the summation 
of wavenumber pairs of long waves meet the Bragg resonant condition (Hasselmann, 
1966): 

2fc/*sin 0 = + (±^ 2 ) (2.13) 
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ojd — iwj + (iu>2) (2.14) 

where ki and k 2 are the wave numbers of the waves. According to perturbation 
mathematical constraints, the second order contribution should be smaller compared 
with the first order Bragg contribution. In practice, some of these higher order 
contributions representing initial manifestation of tilting and multiple scattering 
can become comparable to the first order contribution or larger. Typical relative 
amplitude Doppler spectra were observed by Barrick et al. (1974) using a HF 
frequency Doppler radar (see Figure 2.3). The Doppler spectrum in Figure 2.3 also 
shows a frequency shift of the first order Bragg lines due to advection of the waves 
by currents. 



Figure 2.3: Bragg phenomena in Doppler frequency domain (from Valen¬ 
zuela, 1978). 


Wind may induce surface advection drift, which will add to the Doppler 
shift of the phase speed. Also, rotational (finite amplitude) effects may modify 
the phase speed of short waves. At high winds as the ocean becomes rougher, the 
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higher order contributions become significant, and the Doppler spectra becomes 
broader. The broader Doppler spectra occur when the wind dominant waves are 
longer than the Bragg waves on the ocean. When the Bragg resonant waves are 
generated on longer dominant ocean waves, the mean surface for the Bragg waves 
is no longer flat but undulating. For this typical ocean surface, a composite, or 
two-scale hydrodynamic model is suggested in the literature. 

3. Composite Surface Model 

The observation that Bragg resonant theory by itself does not predict in¬ 
creased radar backscatter in the presence of longer waves led Bass et al. (1968) and 
Wright (1976) to suggest a more tractable two scale-composite surface backscatter 
model. This model represents short Bragg waves riding on longer waves by ap¬ 
proximating the ocean surface as a sum of infinite number of slightly rough patches 
(Figure 2.4). The composite model predicts, therefore, the radar cross-section for a 
single patch as given by (Wright, 1968; Valenzuela, 1978; Brown, 1978): 

<rmi{0i) = Sirk 4 n ghhtfi) + 9wWi) ^(2ka,2k^ sin 6) (2.15) 

where or,- = sin 7 ,-, 7 ; = cos (7 + $)cos 6 , a = sin (7 + $) and 7 ; = cos and 
are the local sea surface slopes (small angles), ghh(0i ) and g vv (0i) are the first order 
scattering coefficients of a local patch, aw is computed with a similar formula, but 
it is much less sensitive to variations in local facet tilts. Thus, the contribution from 
all slightly rough patches over the distribution of slopes p{rj x ,Vv) °f ^ ie dominant 
waves within the radar resolution cell, is given by: 

/ OO f 00 

/ <7////(0i)p(tan $, tan 6)d(tan ty)d(tan 6) (2.16) 

*00 J — 00 

Normally, the backscatter from a slightly rough surface modulated by 
long waves is approximated including both the contribution from a specular point 
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Figure 2.4: Illustration of composite surface approximation (from Stew¬ 
art, 1985). 


model and that predicted by the composite surface model (Hasselmann et al, 1985). 
The tilted facets significantly modify the radar cross-section, especially for the hor¬ 
izontally polarized radiation with a depolarized contribution (Stewart, 1985). 

Towards grazing incidence the polarization ratio decreases as the sea ap¬ 
pears rougher, a phenomenon which illuminated Wright to develop the composite 
model. The cross-section polarization ratio as a function of incidence angle and 
roughness parameter is illustrated in Figure 2.5. At grazing incidence a modified 
model including a combination of composite surface and wedge scattering was sug¬ 
gested by Kalmykov and Pustovoytenko (197C) to explain Doppler spectra from the 
ocean. This combined backscatter model predicts larger backscatter and greater 
Doppler shift for horizontal polarization (at grazing incidence), rather than for ver¬ 
tical polarization. 

4. Application of Composite Model to SAR 

For purposes of interpreting coherent imaging observations, the actual 
scattering from a moving surface like the ocean can be represented by a superposi- 




Figure 2.5: Cross-section ratio due to polarization (from Valenzuela, 
1978). 


lion of the specular and composite theories as described above (however, it is not 
uncommon to include specular scattering in the definition of composite theories), 
plus complex motion effects from a non-rigid distributed surface. The interpretation 
of images from the ocean surface, is still an area of active research. Elachi (1977); 
Hasselmann et al. (1985), and Alpers and Bruening (1986) have given some insight 
to SAR imaging models of the ocean surface. Some of their considerations will be 
described in Chapter III. D. 
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III. PRINCIPAL MECHANISMS OF SAR 
AND INSAR IMAGERY RELEVANT TO THE 

OCEAN SURFACE 

A. INTRODUCTION 

In this chapter the fundamental principles of SAR imaging are briefly de¬ 
scribed first, since the suggested interferometric technique is based on conventional 
SAR. A synthetic aperture radar is an imaging sensor, with a high (two dimensional) 
spatial resolution in azimuth (along track) and in range (across track) directions 
(Ilarger, 1970; Ilovanessian, 1979). High azimuthal resolution is achieved by trans¬ 
lating a coherent radar on-board a moving platform, storing the amplitude and phase 
histories returned from the target, and subsequently processing the stored data to 
remove known quadratic phase errors across the antenna beam. Improved range 
resolution is achieved by using spread spectra techniques, such as pulse compression 
or bi-phase modulation of the radar carrier frequency (Skolnik, 1980). 

The range scanning is at the speed of light in air, but the azimuthal scanning 
is at the speed of the carrier platform (Hasselmann et al., 1985). For processing 
purposes this fact leads to substantial decoupling of the range and azimuth process¬ 
ing. The azimuth channel is related directly to the dynamic distortions associated 
with imaging of ocean surface features, because the time scales of the system and 
the imaged wave fields are of the same order. Inherently, the range channel does not 
have such problems. In this chapter, some fundamental parameters characterizing 
SAR and INSAR operation are described. 
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B. FUNDAMENTAL PARAMETERS 


1. Azimuth Resolution 


At any given time and range, the contributions to the radar image from 


each stationary scattering point at that lange can be identified uniquely by its 


Doppler shift, caused by the sensor motion. This allows the summation of different 


contributions from different scattering points to produce a phase history in the 


aperture synthesis. 


An antenna with a real aperture D has an azimuthal coverage of (see 


Figure 3.1) 



R 


x 


ACROSS 

THE 

TRACK 


ILLUMINATED 
PATCH OP 
SURFACE 


Figure 3.1: SAR geometry. 



(3.1) 


where R is the radar target slant range, and A is the radar wavelength. A R x is the 


real azimuth resolution with antenna of size D. The length of aperture which can 


be synthesized depends on the distance the platform travels while a point scatterer 


is still illuminated. The synthetic aperture for the geometry in Figure 3.1 is A R x , 
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and its maximum azimuth resolution for a stationary rigid point target is: 


A R D . . 

paz0 = 2A R~ x = 1 (3,2) 

the factor of 2 enters because the effective synthetic aperture “beam-width” is half 
of a real aperture. From Equation (3.2), the synthetic aperture resolution is inde¬ 
pendent of the wavelength and range of the radar, and increases with a decrease 
in the physical size of the antenna (Hovanessian, 1979; Robinson, 1984; Stewart, 
1985). 

2. Phase History (of stationary point target) 

Since the resolution in azimuth direction is based on the local Doppler 
shift due to aircraft motion, which in turn is the rate of change of phase, the phase 
history is defined, which is the principle of synthetic aperture imaging radar. Con¬ 
sider a moving SAR with velocity V as illustrated in Figure 3.2. By the time t, the 
radar moved a distance L = Vi, and the range to a stationary point target changed 
from R\ to Ro to R x . Since the phase of the received signal is proportional to the 
distance traveled, the phase of this signal will vary with time. The range R can be 
written as: 

R = [RZ + (Vt ) 7 } l/2 (3.3) 


where R$ is target range of closest approach. By a Maclaurin series expansion 
around t = 0 Equation (3.3) can be rewritten as: 




(3.4) 


Neglecting higher order terms, the phase shift of a stationary point target is (Hov¬ 
anessian, 1979): 


A $(i) = 


2 tt (Vtf 
A Rq 


(3.5) 


19 




Figure 3.2: Top view of range to a point target as a function of aircraft 
velocity (from Hovanessian, 1979). 


Equation (3.5) states that the Doppler shift from a stationary target is a linear 
function of time (chirp). Two additional fundamental parameters which determine 
SAR performance are the processor coherent integration time T and the Doppler 
bandwidth D. 

3. Coherent Integration Time 

The upper limit of the theoretical coherent integration time is deter¬ 
mined by the real aperture azimuth resolution or more precisely, by the real antenna 
beamwidth: 


A R x _ pR 
V ' V 


(3.6) 


where ft is the real aperture beamwidth. 

4. Doppler Bandwidth 

Consider the scenario depicted in Figure (3.3). The Doppler shift of a 
stationary point target B is: 


Ib - 


2V cos 0 

X 


(3.7) 
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Figure 3.3: SAR geometry with respect to azimuth resolution element 
(after Hovanessian, 1979). 


and the Doppler shift for stationary point target A is: 


2Vcos (0 + 0) 


Using small angle approximations for the beamwidth /?, the available Doppler band¬ 
width is: 

21 / 

(3.8) 

5. Time Bandwidth Product 

Combining Equations (3.6) and (3.8) the time bandwidth product is: 



(3.9) 


Equation (3.9) implies that TB > 1, a property which allows significant simplifica¬ 
tions in analysis where the assumption of stationary phase is allowed (Kerr, 1951). 



Figure 3.4: Azimuth image shift. 


C. EFFECT OF POINT TARGET MOTION 

Raney (1971) and Alpers and Bruening (1986), among others, give a clear 
description of the effect of point target motion on SAR image. Following Alpers 
and Bruening (1986), consider a point target with radial velocity U T in the cross 
track (range) direction as depicted in Figure 3.5. The SAR carrier is moving at V 
velocity and R is the distance between the SAR antenna and the point target. 

A SAR which is tuned for stationary targets, locates targets in azimuth at 
zero Doppler frequency. Hence, a stationary reasonable target will be located at 
right angles to the carrier path. However, for the moving target (Figure 3.4), the 
zero Doppler frequency position, will be shifted in azimuth at an angle a = 
relative to a position of a stationary target. Thus the azimuth image shift of a 
target with radial velocity is: 

Ax *s —^ U t (3.10) 

Evidence of this shift is found in SAR imagery for moving rigid targets with velocity 
component in the range direction, like ships or trains. Their images are shifted in 
azimuth direction compared to their wakes or tracks. Keeping in mind fundamental 
principles of SAR, the imaging of the moving ocean surface is discussed next. 











Figure 3.5: Radar hydrodynamic and straining modulation transfer func¬ 
tion and SAR velocity bunching azimuthal transfer function as a function 
of azimuth angle (from Alpers et al., 1981). 

D. APPLICATIONS TO IMAGING OF MOVING 
OCEAN SURFACE 

As described in section B of this chapter, microwave scattering from the ocean 
surface for conventional radars is modeled to an acceptable extent. However, when 
dealing with imaging radars, two new parameters improving the radar performance 
and in a way complicating the backscattering process are introduced: 


• High spatial resolution (order of tens of meters), 
t Doppler information for image generation. 
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The full implication of these two additional factors on the generation of ocean sur¬ 
face imagery is still in active research (Elachi and Brown, 1977; Kasilingham and 
Shemdin, 1988). 

For a stationary surface, a SAR image is a spatial two-dimensional presenta¬ 
tion of the variability of the local coherent backscatter cross section of the surface 
(Elachi and Brown, 1977). For the ocean surface which is nonstationary, SAR image 
is a map of the covariance of complex reflectivity (Hasselmann et al., 1985; Kasil¬ 
ingham and Shemdin, 1988). Due to the difficulties in interpreting SAR imagery 
from the ocean surface, two major resolvable long wave imaging mechanism are 
still recognized, though they represent conflicting view points. One theory is the 
velocity bunching approach, which is based on the coherent particle orbital velocity 
of the dominant wavefield (Elachi and Brown, 1977; Alpers and Rufenach, 1979; 
Swift and Wilson, 1979; Keller and Wright, 1983; Hasselmann et al., 1985; Alpers 
and Bruening, 1986). The other theory is the distributed surface theory which pre¬ 
dicts that the SAR image is composed of amplitude modulation due to straining of 
short backscattering waves and a phase modulation due to tilting by long waves. 
It assumes that the ocean surface is correlated over distances larger than the SAR 
intrinsic resolution. (Jain, 1978; Harger, 1986). 

1. Velocity Bunching Imaging Mechanism 

The velocity bunching theory is based on the fundamental assumption 
that the scattering surface can be divided into individual and uncorrelated scattering 
elements whose dimensions are small compared to the SAR resolution cell (~ 10—20 
m) and large compared to the radar wavelength (5-7 times). Each scattering element 
(facet) is associated with the local orbital velocity and acceleration of the long 
underlying waves. This facet motion gives rise to alternating concentrating and 
spreading of the position of the scattering facets causing amplitude modulation 
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phase locked to the dominant waves. Following Vachon et al. (1988), the radial 
line of sight velocities of a monochromatic long wave associated with a scattering 
element is given by 

U(x,<) = a»/l(sin 2 7 sin 2 <f> + cos 2 7 ) , / 2 sin(k • x — lot + tg~ x {tg^ cos <f>)) (3.11) 

where <f> is the angle between the aircraft velocity vector and the wavenumber vector 
k on the ocean surface plane, and A is the wave amplitude. For swell with a 
period much longer than the SAR integration time, a scattering element originally 
at azimuthal coordinate x will be shifted to (Raney, 1988) 

x' = x + d sin(&cos <f>x) (3.12) 

A constant phase term was not included in Equation (3.12). c' is the 
fractional amplitude of the azimuthal image shift of a scattering element and governs 
the mapping character, c' is defined as (Alpers, 1983) 

c ' = ^A(sin 2 7sin 2 <j> + cos 2 7)^ 2 (3.13) 

Under certain situations the processed SAR image can appear with 
bright intensity patterns which have twice the spatial period of the actual long 
wave on the ocean, which can confuse the interpretation of ocean wave analysis. 
This critical situat ion is quantified by the parameter c and becomes critical for 

c' T = L/2ircos<f> (3-14) 

for c' > L/2ncos<j> the scattering element redistribution becomes many-to-one and 
results in two bright high reflectivity bands associated with one spatial period of 
the actual ocean wave. For c' < L/2ncos<f>, the scattering element redistribution 
is one to one and results in a single unambiguous high reflectivity band related to 
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each spatial period of the actual ocean wave. The parameter c' T is defined by Alpers 
and Rufenach (1979) as a measure of transition between linear and nonlinear SAR 
imaging (see Figure 3.4). In terms of the significant ocean wave height, the critical 
imaging occurs when 


//A 3/2 _ V _ 

\ 7 r / g l ! 2 R cos <£(sin 2 7 sin 2 (f> + cos 2 7) 1 / 2 


(3.15) 


In the above velocity bunching model the ocean surface was assumed commonly 
(Vachon et al., 1988) to be invariant in its structure and scatter dynamics during 
the SAR image formation period (coherent integration time). The integration time is 
considered to be much less than the dominant imaged, wave period. This might not 
be completely true, especially for airborne SAR. Kasilingam and Shemdin (1988) 
and others suggest that the noncoherent scene motion (phase velocity) over the 
image formation period should be explicitly considered in the velocity bunching 
model. 


Additional motion effects of the scattering elements that are sensed co¬ 
herently by the SAR radar is the dominant wave acceleration that corresponds to 
the orbital motion and the scattering element decorrelation time scale. Together 
with the translation of wave patterns during the integration time which enters the 
image noncoherently, these additional ocean surface dynamics blurs the image. The 
scene coherence time efTect is discussed in detail in Chapter IV. The acceleration 
defocus is discussed in Chapter VIII. 

2. Noncoherent Motion Effects 

The translation of mean wave patterns of reflectivity during the course 
of integration time enters noncoherently into the image in two ways: through distor¬ 
tion in wave length and direction known as scanning distortion and through image 


26 







smearing (defocus), also known as look misregistration for a multilook imaging sys¬ 
tem. 

For an monochromatic imaged swell, the noncoherent smearing can be 
corrected by adjusting the processor focus perturbation to one half the projected 
wave phase velocity. This operation will enhance the SAR ocean wave image (Jain 
and Shemdin, 1983). 

Scanning distortion is associated with imaging a moving scene from a 
moving sensor. The along track motion of the scattering cross-section pattern mov¬ 
ing at the long wave phase velocity distorts the azimuthal imaged wavenumber. 
This noncoherent distortion is not related to the coherent phase history of individ¬ 
ual scattering elements. If the dominant imaged wave is propagating only in the 
cross-track direction, only a distortion in wave direction will be practically mani¬ 
fested. The degree of distortion in either wave orientation is proportional to the 
ratio between the wave phase velocity and platform velocity C/V . The distorted 
azimuthal wavenumber is taken to be 

k' az = k az ±k^ (3.16) 

The estimation of scanning distortion from two opposite flights provides 
an opportunity to resolve the ambiguity of 180° in propagation direction that is 
inherent to remote sensors. 

3. SAR Two-Scale Backscattering Model 

Hasselmann et al. (1985^ gives a rather complete theory for Bragg scat¬ 
tering using the two scale approach in a review paper. They distinguish between 
two types of two-scale models: 

• The electromagnetic-hydrodynamic two-scale model. 
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• SAR two-scale model. 


The basic difference between these models is the scale of separation between the two 
wave regimes. In both models, the wave field is treated statistically for wave numbers 
larger than the scale of separation and deterministically for lower wave numbers. In 
the electromagnetic-hydrodynamic two-scale model, the scale of separation, is 
defined by the expansion requirements of electromagnetic hydrodynamic theory. The 
scale of separation has to be chosen one order of magnitude smaller than the Bragg 
backscattering wavenumber, h. On the other hand in the SAR two-scale model, 
the scale of separation, k aar , is determined by the inverse of SAR resolution scale. 
The SAR two-scale imaging model is based on the physics of the more fundamental 
two-scale electromagnetic-hydrodynamic model. Typically the fc em /, is one order of 
magnitude larger than k 3ar (see Chapter IV). 

4. SAR Mapping Transfer Function 

In addition to the Doppler induced modulations (velocity bunching) typ¬ 
ical for coherent imaging radars, real radar cross-section modulation tilt and hydro- 
dynamic modulation are commonly used to describe ocean surface imaging. These 
imaging mechanisms govern the imaging process for waves propagating close to range 
(cross-track) direction. In this orientation the velocity bunching is nil (see Figure 
3.4). For moderate sea states the mapping of the wavefield might be described by a 
linear transfer function 


SAR = j^Hydr y^bunching 


(3.17) 


Usually this situation does not occur, since the mapping is non-linear and the super¬ 
position principle does not hold, i.e., the image spectrum is distorted from the actual 
wavefield. Electromagnetic interaction , or tilt modulation of backscatter radar cross 
section, is caused by the actual variation of local incidence angle. Hydrodynamic 
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modulation of backscatter radar cross section is caused by the interaction of the 
Bragg waves with the long waves on which they are superimposed. The long wave 
orbital velocity strains the surface producing areas of current convergence and di¬ 
vergence, resulting in hydrodynamic modulation of the amplitude and wavenumber 
of the Bragg waves. The surface displacement distorts the air flow over the sur¬ 
face, causing short waves to grow at varying rates (modulated wind stress) for a 
given wind speed, depending on their relative position on the long wave (private 
communication, Valenzuela, 1989). 

E. INTERFEROMETRIC SAR (INSAR) 

1. Principle of Imaging of Ocean Surface Velocity Fields 

Interferometric synthetic aperture radar (INSAR) is a relatively new 
method for measuring ocean surface wave fields, spectra, and currents. Unlike 
conventional SAR, INSAR has the potential for measuring unambiguously virtual 
surface velocity fields. INSAR combines large coverage with high nominal spatial 
resolution (~ 12 m in this work), good sensitivity and appears suitable for both 
airborne and probably satellite implementation. The airborne INSAR technique 
has been developed at Jet Propulsion Laboratories (JPL) by Goldstein and Zebker 
(1987), and was installed first on a NASA CV990 aircraft and recently on a DC-8 
for evaluation experiments. 

Following Goldstein and Zebker (1987), the basic technique consists of 
aligning two antennas along the radar flight path. Pulses are transmitted from the 
rear antenna and are received simultaneously by both antennas. The resulting sig¬ 
nals are converted, separately, into two complex images in the conventional manner 
for SAR. Any line of sight motion of the identical imaged surface (resolution cell) 
observed by both antennas is translated into phase difference when these resolution 
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cells are combined interferometrically. This phase difference, cr, is proportional to 
the virtual measured radian Doppler shift, Athat occurs in the time it takes for 
the rear antenna to move to a place of symmetry with respect to the front antenna 
(see Figure 4.2), where 

u D = 2k • U (3.18) 

and 

a = WjD At = y y U r (3.19) 

where B is the baseline separation between antennas, V is the aircraft velocity, U r is 
the radial c. mponent of ocean surface motion and A is the radar wavelength. Equa¬ 
tion (3.19) does not consider speckle contribution from the random ocean surface 
and/or from the measurement and/or processing system errors. These limitations 
are reviewed briefly in the next section and further discussed in the next chapter. 

The resulting interferogram is a complex image with phase which is the 
difference of the two original SAR image phases and magnitude which is the product 
of the two original magnitudes. Practically, the phase and amplitude images are 
processed into two separate images. For more details see Chapter IV. B. 

2. INSAR Limitations 

The limitations of INSAR for spatial imaging of ocean surface velocity 
structure may be categorized into two types: 1) system dependent 2) the degree 
of coherency of the imaged ocean scene within the inherent time measurement of 
INSAR. 

a. Platform Stability 

In order to compare the phases of a rezel (resolution cell) which is 
observed at two different times, the viewing aspects must be nearly the same. This 
might not be fulfilled due to lack of accurate knowledge of relative antenna location, 
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which in turn will result in non-zero velocities and acceleration in the radial direction 
(antenna pattern disimilarity). For Equation (3.19) to extract accurately the radial 
component of ocean surface motion, it is assumed that the two antennas are collinear 
with the radar motion. Aircraft pitch and yaw will add phase difference between 
corresponding pixels, since one antenna will be closer in range to the observed 
resolution cell. The radar cannot distinguish between a slight sideways drift of 
its carrier and the radial motion of the ocean surface. This contamination can be 
compensated using readouts of aircraft inertial navigation system (INS), or to a 
reasonable degree, by inclusion of land in the scene to give absolute reference, or 
a combination of the two. The INS provides aircraft pitch, yaw, pitch and yaw 
rates, horizontal and vertical velocity, and other flight parameters. With these 
data, the expected error phase shifts from a stationary scene can be computed and 
compensated accordingly. A yaw error of 0.01 degrees will yield an error in current 
estimates of 4 cm s _1 (Goldstein and Zebker, 1987). 

b. Backscattering Surface Roughness Coherence 

The interferometric technique requires comparison of two samples 
of a moving surface element (resolution cell) in a time interval A t = B/2V required 
for the rear antenna to move (the mean aircraft speed is 220 m/s) half the distance 
between the antennas B/2 = 9.9 m. It is assumed that in this time increment, A t 
= 0.041 s, the surface wave field as a whole in the resolution cell (~ 6 -r 8x 12 m) 
does not decay or grow, i.e., the scattering patch is nearly stationary. A number of 
hydrodynamic processes afTect the lifetime of Bragg resonant waves. Among these 
are viscous dissipation and their response (relaxation) to external forcing, such as 
wind, currents, resonant nonlinear interactions, straining by orbital velocity of long 
waves, and turbulence. The viscous half life time decay of a 10 cm wave is 63 s. But 
the lifetime of wind waves are controlled by relaxation processes, which are inversely 
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proportional to the wind growth rate. Therefore, as the wind speed increases, the 
relaxation time decreases. For 10 cm waves, the relaxation time is about two seconds 
for 6 m/s winds (Valenzuela and Wright, 1976). Hence, for moderate sea conditions 
the ripples should remain nearly constant within the time interval used to measure 
the phase difference. However, the spread of orbital velocities in the radial direction 
within the resolution cell will also decorrelate the backscatter complex signal causing 
an ambiguity of the phase as the ocean wave height changes. An initial attempt to 
study the limitations of the INSAR technique due to finite scene coherence time is 
made in the next chapter. 
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IV. THEORETICAL ANALYSIS OF INSAR IN 
IMAGING OCEAN SURFACE DYNAMICS 


A. INTRODUCTION 

In Chapter III. D, the known physical causes for the ocean surface decorre¬ 
lation were described. The chapter ended with a recommendation that the velocity 
spread of the scattering elements due to local orbital velocity of the subresolution 
scale waves should be further analyzed. In this section for the first time, a prelim¬ 
inary theoretical model for an airborne INSAR imaging the ocean surface is given. 
This model is based on the accepted SAR theory of imaging the moving, nonrigid, 
distributed ocean surface (e.g., Alpers and Rufenach, 1979; Swift and Wilson, 1979; 
Valenzuela, 1980; Raney, 1981; Jain, 3 081; Rotheram, 1983; Plant and Keller, 1983; 
Tucker, 1985; Hasselmann et al., 1985; Monaldo and Lyzenga, 1986). It is shown 
that the limited ocean surface temporal coherency due to the velocity spread of 
SAR subresolution waves might limit the capability of INSAR in providing virtually 
instantaneous high resolution snapshots of near ocean surface velocity patterns. 

In the analysis given below the SAR two-scale model is assumed implicitly. 
This model was introduced first by Tucker (1985) who modified the earlier two-scale 
hydrodynamic model (Wright, 1968, Valenzuela, 1978). The model divides the ocean 
surface waves into two scales, a subresolution short scale with wave numbers greater 
than the SAR azimuthal cutoff wavenumber which are treated statistically, and a 
long scale with wave numbers shorter than this wavenumber which are treated deter¬ 
ministically (see Chapter III. D). It is found that the motions due to subresolution 
waves produce an azimuthal smearing which can be represented by a Gaussian low 
pass filter acting on the azimuthal component of wavenumber in the image (Tucker, 
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1985). This degradation in azimuth resolution (smearing) is inversely proportional 
to the scene coherence time and determines the nominal azimuth resolution for SAR 
and INSAR. The scene coherence time for L-band radar will be estimated based on 
the Pierson-Moskowitz ocean spectrum (350-600 msec) and the backscatter reflec¬ 
tivity of simultaneous SAR and INSAR observations from the Marina experiment 
on September 8, 1989 (~ 110 msec). 

A major accepted simplification of SAR imaging theory, both for mean image 
and speckle statistics, follows from the hypothesis that the complex reflectivities at 
different points on the sea surface are uncorrelated. The spatial separation of the 
different scattering elements (facets) should be large in comparison with the Bragg 
wavelength, but still small in comparison with the SAR resolution scale. The surface 
reflectivity may therefore be regarded as spatially white (after Hasselmann et al., 
1985). This assumption will be used in the dirivation of the preliminary model of 
the ocean surface cross-covariance of complex reflectivity for INSAR. 

B. INSAR IMAGING OF THE OCEAN 

The complex radar return signal at time t is given by (Kasilingam and 
Shemdin, 1988): 

S(t)=f f r(x,t)f(t-2 Rq/c)A(x-V t,y)dxdy (4.1) 

with x = x(x,y), see Figure 4.1. 

where r(x, t) is the time variant ocean surface Bragg and specular reflectivity, f(t — 
2R 0 /c ) is the radar pulse function and A(x — Vt,y) is the antenna pattern. We shall 
ignore the compression effect in the pulse function since the time variation of the 
returned signal within the pulse duration (few microseconds) is negligible. The range 
dependence of the pulse function introduces a quadratic phase term. For simplicity 
assume that the pulse rate is high enough to fulfill the sampling criteria in azimuth 
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and a continuous signal is transmitted. Therefore, only a quadratic phase term will 
be introduced by the pulse function. An illustration of an airborne SAR imaging 
the ocean surface is depicted in Figure 4.1 where k is the radar wave number in 
the viewing direction, V is the platform velocity and Ro is the initial aircraft target 
slant range (see Figure 4.1). 



Figure 4.1: An artistic illustration of SAR geometry in imaging the ocean. 

The recorded received signals are processed by convolving the records with a 
matched filter, resulting in the SAR image. For a focused system (see Chapter III) 


35 




with no focusing adjustments for wave pattern motions, the matched filter impulse 
response is taken to be 

h(x - Vt ) = exp [tfc(x — Vt) 2 /Ro] (4.2) 

The complex conventional SAR image of a pixel from a single antenna is 
therefore obtained by convolving S(t ) given by Equation (4.1) with the matched 
filter (Equation (4.2)): 

/ ■fee /4*00 /•-+* 00 

S{t — r)h(T)dr— / / r(x,t - r) 

•oo J—oo J-oo 

•exp [— ik(x -V(t- r)) 2 /f?o] A(x — V(t — r)) 

•exp [ijfc(x - Vr) 2 /f?o] dxdr (4.3) 


where the first exponential term on the right hand side of Equation (4.3) is the 
quadratic phase history for a stationary target. Similarly, the corresponding complex 
pixel image for the second antenna taken at a slightly different time A t = Bj 2V, 
which is the time required for the rear antenna to move to a place of symmetry with 
respect to the forward antenna, is given by: 


/ +oo 

S{t + B/2V - T , )h{r , )dT t 

•00 

/ +00 f+OO 

/ r(x ,t + B/2V — r) 

-oo J—oo 


•exp [ik{x - V(t + B/2V - t')) 2 /R]A(x - V{t + B/2V - r ')) 
•exp \ik(x — Vt') 2 /R^ dx!dr' (4.4) 


Equations (4.3), and (4.4) represent complex SAR images of the ocean surface 
reflectivity of the same pixel. The goal is to extract reliably the mean Doppler shift 
caused by the first order ocean surface motion. This phase information is included 
entirely in the cross-covariance of the complex surface radar reflectivity (rjr^) (see 
Equation (4.6) below). From a statistical analysis point of view, the averages of 
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rir 2 (denoted by brackets) are defined for an infinite ensemble of different surfaces. 
Practically, the averages are estimated over a finite surface region for a finite time. 
By taking the cross-covariance of the corresponding complex images from the two 
antennas, the phase difference, which is proportional to the virtual surface velocities, 
is obtained. The interferometric operation is therefore given by: 


Q( X ) = (tj(< = x/V ) Q(< + B/2V)) (4.5) 


Substituting Equations (4.3), (4.4) into (4.5) yields: 

r+oo /+oo f+oo /+oo 


/ +oo r+oo /+oo f+OO 

III (ri(x,<-r)r;(x # ,< + B/2V - t')) 

*00 J—oo J -CO J—oo 

•exp [— ik(x — V(t — t)) 2 /Ro] A(x — V{t — r)) exp [ik(x - Vr) 2 /^] 
•exp [tfc(x' - V(t + B/2V - T'WlR^Atf - V(t + B/2V - r')) 

•exp \ik(x' — Vt ') 2 jRo]dxdx!dTdr' (4.6) 


The first term on the right hand side of Equation (4.6) is the covariance of 
complex reflectivity of the ocean surface, which is defined by the second moment 
statistics of the ocean surface backscattered signal. As noted in the introduction 
to this section, the averaged SAR mean image is directly related to these statistics, 
namely, the covariance of complex reflectivity of the ocean surface. Figure 4.2 
illustrates the geometry of INSAR. 

The extraction of the phase information from the backscattered moving sur¬ 
face by using Equation (4.5) is the essence of INSAR in imaging the ocean. This is 
fundamentally different from the SAR image intensity, given by: 


(i*i P) = 


(4.7) 


Equation (4.7) clearly illustrates the loss of the phase information in the SAR image, 
which is retained in Equations (4.5), (4.6) for INSAR. 
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Figure 4.2: Geometry of INSAR: a) top view, b) perspective view shows 
the change in distance of the resolution element due to wave motion in 
the time it takes for the rear antenna to more to a symmetric position 
with respect to the front antenna. 

The SAR System Transfer Function (STF) is not considered at the present 
stage. It is assumed that the SAR is a perfect instrument that maps the covariance 
of complex reflectivity into a two-dimensional complex image. The interest here is to 
investigate the feasibility of INSAR in imaging the moving ocean surface. Therefore, 
the interferometric operation described by Equation (4.5) can be replaced by the 
cross-covariance of complex reflectivity of an imaged ocean scene sensed by the two 
INSAR antennas with a time difference At = B/2V. 

C. DERIVATION OF THE CROSS COVARIANCE OF 
THE COMPLEX RADAR REFLECTIVITY FOR 
INSAR 

The complex radar reflectivity for an individual scattering element moving 
with velocity U(x,t) is given by: 

r(x, t) = r 0 (x, t) exp|2ik • U(x,t) tj (4.8) 
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where ro(x, t) is the instantaneous complex reflectivity of an ocean surface scattering 
element within the SAR resolution cell. As noted earlier, these scattering elements 
(facets) are assumed to be large relatively to the resonant backscattering Bragg 
waves (5-10 times), but much smaller than the SAR resolution cell (rezel). 

The cross-covariance of complex reflectivity of a time dependent backscat¬ 
tering surface f::om corresponding rezels as imaged by the two antennas is given 
by: 


(ri(x, t ) rj(x', t')) = (r 0 i(x,i)exp [2tk • U(x, t) <] r^x', t') 


•exp 


-2ik • U(x', t ') t'| )S(x - x') ctq(x, t) (4.9) 


where t' — t + At = t + B/2V and S(x — x') represents the white spatial distribution 
of the scattering elements within the SAR resolution cell (Hasselmann et al., 1985). 
<7 0 (x,t) is the variance of r and represents the modulation of the covariance of 
short wave reflectivity by the slow time varying long wave tilt and hydrodynamic 
modulation. cr 0 (x,t) represents the radar cross-section as imaged by conventional 
SAR and real aperture radars. 

Ir Equation (4.9), the t me lag At represents only the constant time interval 
required for INSAR measurement. The additional time difference of r — t' (see 
Equation (4.6)), which is of the order of the SAR integration time, is not considered 
here explicitly (except when determining the effective azimuth resolution). This 
time difference is not relevant for INSAR operation once the two corresponding 
complex SAR images are formed. For more detail on the covariance function for 
SAR, see for instance Hasselmann et al. (1985). The discussion here (Equation 
(4.9)) is based on analyzing the second order statistics of the fluctuating reflectivity 
represented by the ensemble term in Equation (4.9) on the time scale relevant for 
INSAR (At = B/2V). 
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The SAR two-scale model is invoked resulting in the ocean surface velocities 
decomposed as 

U(x, t) = U 0 (x, t) + U'(x, t) (4.10) 

where U 0 (x, t) is the long scale velocity and U'(x, t) is the subresolution short scale 
velocity. Uo(x,t) can be extracted by INSAR if the temporal variability of the 
backscattering surface, defined by the radar scene coherence time, is small relative 
to the time difference At = B/2V required to obtain the phase image of INSAR 
which is directly proportional to the mean radial surface velocity. 

The limits of the temporal decorrelation of the ocean surface on INSAR in 
obtaining the surface velocity U 0 (x, f), are now addressed. The variation of the 
ocean surface velocity U(x,f) in the time At = B/2V can be described by using 
the Taylor expansion 

U(x,< + B/2V) = U(x,<) + B/2V (4.11) 

Substituting Equations (4.10) and (4.11) into Equation (4.9) and factoring out the 
deterministic long scale contributions yields 


(ri(x,<) rZ(x,t + B/2V)) 


exp —2ik-U 0 (x,f) R/2Vj 
• exp f—2ik 


^Ho(x t )± 

dt {X ' t} 2V 


roi{x,t) r^ 2 (x,t + B/2V) 


•exp 


-2i^k(u'(x,«) + 



Equation (4.12) represents the fundamental basis of the present study in evaluating 
the feasibility of INSAR in imaging the ocean surface. It combines the relevant 
physics of the backscattering surface as a function of INSAR physical configuration. 
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The second exponential term on the right hand side of Equation (4.12) can 
be neglected if 

^^( (+B /2V)<l (4.13) 

where t is of the order of the SAR integration time. Since t > B/2V and the 
absolute value of the acceleration of the dominant wave is giver by d\Jo/dt = wUo, 
inequality (4.13) can be rewritten as: 

wl« 1 (4.14) 

or 

T>2nt (4.15) 

where T is the period of the imaged dominant wave. Equation (4.15) specifies the 

limiting condition for the integration time when imaging dominant waves by SAR 
and INSAR. The effective processing integration time for mean incidence angles 
(30°) is almost 1.7 sec. In practice, the scene coherence for a given radar wavelength 
dictates the effective integration time, which is substantially less than 1 sec (e.g., 
Tucker, 1985). This crucial issue is discussed in the next two sections. 

The first term on the right hand side of Equation (4.12) is the required phase 
shift to be determined by INSAR. It represents the first order ocean surface radial 
velocities. U 0 is a vector sum of 

U 0 = U cp + U ob + U cr + ^ar (4.16) 

The individual contributions are identified as: 

• U cp phase velocity of the backscattering resonant Bragg waves; 

• U ob mean orbital velocities of waves substantially longer than the character¬ 
istic dimension of the resolution cell. 
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• Uch material velocity of any near surface currents underlying the backscat- 
tering Bragg waves. 

* Uytfl aircraft attitude error velocity due to yaw and pitch bias of the INSAR 

antennas (see Chapter III. B). 

♦ 

Practically, the ensemble average factor in Equation (4.12) expresses the ad¬ 
ditional fading in the cross covariance of complex reflectivity due to the slight time 
difference A t when the same scene is imaged separately with a small time delay 
by the two INSAR antennas. This phenomena is defined here as c(x, t, At). These 
modifications in Equations (4.13), (4.15) substituted into Equation (4.12) yield: 

(ri(x,t) rj(x,t + B/2V)) = exp[-2ik-U 0 (x,<) ■£?/2V’]c(x,t, B/2V)ao(x,t) (4.17) 

Both are shown to be quantified as a function of the finite scene coherence 
time. c(x,t, At) can be rewritten as (see Equation (4.12) (4.17)): 

c(x,« > B/2V) = expf-2<^k.(u'(x,i) + ^^ (< + ^))j </ ( x ’ i ’ B / 2V '> 

(4.18) 

with 

c'(x, t, B/2V) = (r 01 (x, t) (x, t + B/2V)) (4.19) 

where d(x,t,B/2V) is the temporal cross-correlation function value for the INSAR 
time difference At = B/2V. The phase term in Equation (4.17) can be rewritten in 
terms of radian Doppler frequency 

lo, s = 2k • U'(x, t) (4.20) 

where u 3S is the perturbation in radian Doppler shift associated with the velocity 
spread of modulating the subresolution scale wave scattering elements. Substituting 
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Equation (4.20) into Equation (4.18) yields: 


c(x,t,B/2V ) = exp 




(4.21) 


where du> t »/dt ( B/2V ) 2 , which is 0(1O~ 2 ), is neglected. 

In contrast with SAR image intensity which ignores phase information (Equa¬ 
tion (4.7)), INSAR imagery represents the virtual phase difference between the two 
corresponding SAR images. Therefore, the present model (Equations (4.14-4.21)) 
takes into account additional complex terms which were not considered explicitly in 
SAR models. We shall further concentrate in evaluating the complex disturbance 
to INSAR in sensing ocean surface velocities introduced by c(x, t, At). This dis¬ 
turbance includes cross section attenuation and additive phase noise, both caused 
by the subresolution intermediate waves modulating randomly the backscattering 
elements (facets). The finite temporal coherence of the ocean scene is of special im¬ 
portance since it defines the maximum coherent radar integration time for a given 
radar wavelength, and therefore the maximum azimuth resolution achievable. In the 
p. esent study, the effect of the finite scene coherence in the time interval ( B/2V ) 
required to obtain the interferometric information based on the corresponding inte¬ 
grated SAR images from the two antennas is discussed. 

When considering the radar scene coherence time, two time scales should be 
looked at: 


• The decorrelation time of an individual facet within the SAR resolution cell. 
This fading in the mean complex reflectivity can be described by a temporal 
decorrelation function. 

• Intrinsic coherence time in which the Bragg waves of a single scattering el¬ 
ement (facet) preserve their structure. It is assumed that this time scale is 
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long compared with the coherent integration time (Hass., 1985, Vachon et al., 
1988). 


The temporal correlation function is the inverse Fourier transform of the 
reflectivity variance Doppler spectrum. Based on theory (e.g., Tucker, 1985) and 
the TOWARD experiment (1984), it can be assumed that this Doppler broadening 
of short wave reflectivity has a Gaussian shape. After removing the mean Doppler 
shift due to long wave and current velocities, the velocity power spectrum of the 
short wave reflectivity from a given pixel is 

P{u) = J ( r 0 i(t ) rj 2 (< + r)) exp[-iu;T]dr = P 0 exp [-u; 2 4 /2(7 2 ], (4.22) 

where a u is half Doppler bandwidth of P(u>). 

The inverse Fourier transform of Equation (4.22) yields the temporal cross¬ 
correlation function of short wave reflectivity at a given point. 

1 r+°° r t 

c'(t) = — / P(u>)exp [iwrjdw (4.23) 

Z 7 T J-oo *■ ■* 


Substituting Equation (4.22) into Equation (4.23) yields 



+oo 

exp 

-oo 


-u; 2 /2^ 


exp 



(4.24) 


where u = u> aa . 

Performing the integration in Equation (4.24) renders 

At) = exp [SS/i] = ^ exp [-^r 2 /2] (4.25) 

The standard deviation of the Doppler broadening of short wave reflec¬ 
tivity represents the inverse of the correlation time of the subresolution scattering 
elements (facets), namely, the finite scene coherence time. 


(4.26) 
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Substituting Equations (4.25) and (4.26) into Equation (4.21) yields 
c(x,*,r) = P { l exp UJ exp [ir (u„ + , 


c ( x >^>r) = Pq exp + ( 4 - 27 ) 

where Pq = <r u ,Po(27r) -1 ^ 2 and r = B/2V. 

The phase term in Equation (4.27) expresses the Doppler spread of the scat¬ 
tering element (facet) bandwidth, which is assumed to be the same for all facets. 


Thus 


l w «# + 


du>., t 


O(Au) = 


Substituting Equation (4.28) into Equation (4.27) yields 


, . v D , r 2 l \i2nr 

c(x,t,T) « P Q exp — r— exp - 

[ 2r c J r c J 


(4.28) 


(4.29) 


Equation (4.29) represents explicitly the previously mentioned complex dis¬ 
turbance to INSAR remotely sensing the ocean surface. Substituting Equation 
(4.29) into Equation (4.17) yields the explicit expression of the cross-covariance of 


complex reflectivity for SAR interferometry. 


(rj(x, t) At)) = cr(x,t) exp - 

L t c J 

•exp — exp[—2ik• Uo(x,t) Atj (4.30) 


•expj^-—j exp[-zm-u 0 tx,tj arj 

where At — B/2V. Equation (4.30) illustrates the sensitivity in measuring the line 
of sight velocity component U 0 of an unsteady rough ocean surface. It is shown 
that the key limiting parameter is the ratio between the inherent interferometric 
time delay (At = B/2V) and the scene coherence time t c which is analyzed next. 

The first phase term in Equation (4.30) represents the additive phase noise to 
the mean surface velocity Uo, which is attempted to be measured by INSAR. This 
noise term is required to be small relative to the mean phase term. This implies 


|2k • U 0 At| > 


I 27tA 1 1 


(4.31) 






which results in 


Tc> 2|U 0 |' 

For the NASA/JPL DC-8 airborne configuration, the range of measurable velocities 
Uo is ±1.25 m/s and for the L-band radar wavelength of A = 24 cm, Equation 
(4.32) indicates that the scene coherence time r c should be at least 0.1 sec for the 
additive phase noise term to be neglected. The known observed estimates of the 
scene coherence time for L-band radar are equal find smaller than the order of 10 2 
msec (Lyzenga and Shuchman, 1983, Kasilingam and Shemdin, 1988). 

Next to be examined is the attenuation of the radar cross-section caused by 
the slightly different time in imaging a given scene by the two physically separated 
INSAR antennas. This attenuation is represented by the second exponential in 
Equation (4.30). The radar cross-section for INSAR can be expressed as: 

</(x, t , At) = <ro(x, t) exp(—Af 2 /2r 2 ) (4.33) 

In order to extract the phase information representing the surface velocity 
(last term on RHS of Equation (4.30)), one needs a sufficient signal <r'(x,f,Af) 
above the noise level. A 3 dB attenuation of the cross-covariance ‘signal’ is allowed, 
ensuring accurate measurements of the backscattering surface velocity. The required 
scene coherence time r c in Equation (4.30) is taken to be 

Te ~ 2V (2 In 2) 1 / 2 ( 4 ‘ 34 ^ 

Once t c is known, the upper limit on At = B/2V is defined, which in turn prescribes 
the necessary INSAR set up. Conversely, for a given INSAR set up the minimum 
allowable scene coherence time is determined. For the DC-8 INSAR configuration 
(B = 19.8 m, platform velocity V « 200 m/s), using criterion (4.34) yields the 
minimum required scene coherence time r c = 20 msec. Combining the phase and 
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amplitude constraints on the radar scene coherence time, the temporal correlation 
of the short wave reflectivity should be substantially greater than 100 msec for the 
given INSAR system. 

Finally, the model can estimate the range of variability of the complex re¬ 
flectivity within the inherent INSAR measurement time, which will still allow the 
extraction of ocean surface resolvable motion with minimum contamination. 

D. ESTIMATION OF SCENE COHERENCE TIME 

The next step in the study is to estimate the actual scene coherence time for 
an L-band imaging radar using two different approaches: 

• Theoretically, assuming the Pierson Moskowitz ocean spectrum. 

• Experimentally, using for the first time simultaneously SAR and INSAR imag¬ 
ing records. 

1. Estimation of Scene Coherence Time Based on 
Pierson-Moskowitz Spectrum 

The finite scene coherence time may be one of the most important factors 
in constraining SAR’s ability to image azimuthally traveling waves (see Section IV. 
B). In this section, the analysis will follow recent attempts to quantify this parameter 
(see e.g., Tucker, 1985b, Alpers and Bruening, 1986, Vachon et al., 1988) using the 
Pierson-Moskowitz ocean spectrum for fully arisen sea. Kasilingam and Shemdin 
(1988) point out the uncertainty in the validity of this approach, which might lead to 
differences by order of magnitudes in estimating the scene coherence time r c . Due to 
the complexity of the temporal ocean surface decorrelation r c , the theory associated 
with the ocean scattering motion effects is still under extensive research. Owing to 
the lack of experimental evidence and more profound theoretical understanding of 
the small scale roughness of the sea surface, this approach still remains acceptable. 
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The classical approach is to relate the azimuth image smear 
(6(x) 2 ) 1 / 2 (which is inversely proportional to the scene coherence time), to the rms 
facet surface velocities (U ' 2 ) 1 / 2 caused by the random material velocities of subres¬ 
olution scale waves (see Equation 4.20). 

{6(x) 2 )'' 2 = (V ,2 )^ 2 R/V (4.35) 

where {U ,2 ) x ^ 2 is the radial component of (U' 2 ) 1 / 2 . 

The shape of the Doppler broadening spectrum associated with the ran¬ 
dom facet velocities U' is assumed to be Gaussian (see Equation (4.22)), resulting 
in a SAR azimuthal imaging transfer function with a low-pass filter character. The 
bandwidth of this filter is inversely proportional to the decorrelation time scale, 
which is interpreted as the scene coherence time. In the wave number domain, the 
width of this filter (transfer function) is inversely proportional to the azimuthal im¬ 
age smear. It is important to note that the loss of resolution due to finite scene 
coherence time is not a defocusing effect, and is therefore not recoverable by focus 
adjustment in the processing or integration time setting (see Chapter III). 

Based on the relation between the scene coherence time and rms surface 
velocity (U ' 2 ) 1 / 2 spread, the subresolution velocity spread is first estimated. (U ' 2 ) 1 / 2 
may be computed by integrating the orbital velocity power spectrum between the 
appropriate wave frequencies 

( u r 2 ) = S„(u) f 2 { 7 ) <Ju>, (4.36) 

JUf\ 

where S v {u)) is the omnidirectional surface elevation frequency power spectrum and 
/( 7 ) is the geometric function 

fil) = ~ / (sin 2 7 sin 2 rf> + cos 2 jy^ 2 d<f> (4.37) 

7T J-1 r/2 
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where <f> is the angle between the wave train and the aircraft attitude and 7 is the 
incidence viewing angle (Figure 4.1). The geometric function /( 7 ) represents the 
projection of the moving scene in the radial viewing direction. /( 7 ) is always smaller 
than 1 . 

A rather crude assumption is next made about the ocean wave spectrum 
in calculating (U' 2 ) 1 / 2 , which is in fact sensitive to local sea conditions. In the 
absence of relevant experimental data, it is assumed that the ocean wave spectrum 
for a fully arisen sea can be described by a Pierson*Moskowitz spectrum (e.g., Alpers 
and Bruening, 1986, Tucker, 1985b). Substituting the Pierson-Moskowitz spectrum 
into Equation (4.36) yields: 

(U' 2 )=jf tu 2 a< 7 2 uT 5 exp [—/%/U w ) 4 ] f 2 {l)du> (4.38) 

where a = 8.1 x 10 -3 is the Phillips constant, /? = 0.74, and U w is the wind speed 
at a height of 19.5 m. 

After integration, 

0 1/2 

ui «f \p> (U w )' 2 ] / 2 ( 7 )|” (4.39) 

with p = (Py^g 2 = 82.78 m 2 sec -2 . 

The lower limit of integration is associated with the shortest resolvable 
ocean wavelength 0 (10 m) and is determined by the Nyquist criterion. The upper 
limit is determined by the fundamental subresolution scattering elements (facet) 
with size of the order of Bragg wavelength (~ 0.2 m) which can be replaced for 
computational convenience by u> 2 —♦ 00 with a small error. 

The lower limit is controversial since the actual SAR (and INSAR) az¬ 
imuth resolution is not known a priori, being a function of the scene coherence 
time. A few choices based on different physical considerations will be discussed 
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below. The lower integration limit u>\ in Equation (4.39) is replaced by a cutoff 
resolvable wavelength L\ using the dispersion relation for water waves, 

= (2irg/Li) lf2 (4.40) 


Since the shortest resolvable wavelength is 0 (10 m), the surface tension effects are 
neglected and only gravity waves are considered. By inserting the integration limits 
into Equation (4.39), it is found 

«’>=(!) gyV'* [✓(£)] m <«*> 


with n' = \ijg2-K = 1.35 m sec -2 . For strong winds ^'Li/U 2 < 1, the radial velocity 
spread from Equation (4.41) simplifies to: 


( u ; 2 ) 1/2 = 



m 


(4.42) 


which is independent of the sea-state. For weak wind speeds \i'L\f U 2 > 1 and 
Equation (4.41) simplifies to: 


( u : 2 ) i/2 = 


a 

4 


1/2' 


Uw Z(7) 


(4.43) 


Surprisingly, the rms velocity spread for weak winds, Equation (4.43), is indepen¬ 
dent of the cutoff wavelength Xj, which is determined in this case mainly by the 
SAR system. This wavelength represents the scale limit of waves which are treated 
deterministically (L > L\) and waves which are treated statistically (L < L\). This 
statistical region is of interest since the unresolved wave (L < L\) modulate ran¬ 
domly the fundamental scattering elements. The temporal variability in the surface 
reflectivity is determined by the scene coherence time. 

Two related choices for the cutoff wavelength L\ have been considered. 
The first choice is based upon the processor coherent integration time for a wave 
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with a period equal to T. This wave and shorter waves will propagate at least 
one cycle over the time of image formation contributing only to the rms velocity. 
Using again the deep water dispersion relation for gravity waves, this cutoff scale is 
(Vachon et al., 1988): 

u = far 1 (4.44) 

The second choice used by some investigators (e.g., Alpers and Bruening, 1986, 
Tucker, 1985a) is L\ = 2 p, twice the matched filter processor resolution. For a mov¬ 
ing ocean surface, the azimuth resolution is degraded and is practically unknown. 
It could be much larger than 2 p, which is the nominal azimuth resolution for an 
unmoving scene. Noting that the mean Doppler broadening due to the velocity of 
the intermediate unresolvable waves is given by: 


_ 2( U;*)»/* 

Jdc — X 


(4.45) 



(4.49) 








TABLE 4.1: Estimates for Scene Coherence Time as a Function of Wind 
Speed and Cutoff Wavelength for NASA/JPL DC-8 L-Band Radar with 

fin) = i 


Equation 

lx 

U w (m/sec) 

r c (msec) 

(4.47 ) 1 

JLT 2 

It 1 

>10 

600 

(4.48) 2 

2 P 

> 10 

610 

(4.51) 3 



350 


1 T = 2 sec 2 p = 3.03 m 3 H, = 0.6 m. 


Since for fully developed wind seas, the wind speed U w at the height of 19.5 m is 
related to the significant waveheight (mks units) by 


U 2 = 47.4#„ 


(4.50) 


Equation (4.49) can be rewritten as 


T '~ hl2 VW,M (4 ' 51) 

The estimated values of the scene coherence time in Table 4.1 using 
Equations (4.47), (4.48), (4.51) and /( 7 ) = 1 are rather long compared to the 
sparse estimates based on observations in the literature (30-100 msec for L-band 
radar). The decrease in f(y) with increasing incidence angle will increase even more 
the scene coherence time r c as computed by the equations used in Table 4.1. The 
dependence of r c on the spreading function /( 7 ) is illustrated in Figure 4.3. This 
dependence was not considered explicitly in other works and its importance rises 
from analysis of observed data (see Chapter VII). 

Variability of several hundred msec in the estimates of coherence time 
scales is noted in Table 4.1. The first two estimates are computed using the DC -8 
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Figure 4.3: The effect of spreading function /( 7 ) on the scene coherence 
time r c , showing the increase of the scene coherence time with incidence 
angle. 

SAR parameters and are independent of sea state, a situation which seems to be 
unlikely. The third estimate seems to be more sensitive to the real ocean situation, 
namely, dependence of scene coherence time on sea state and hence wind speed. It 
will be shown through the experimental approach adopted in the next section that 
the shortest coherence time estimate of the present method is rather longer than 
that which is observed. This might indicate the inadequacy of the present approach 
to estimate the scene coherence time, mainly because of a lack of comprehensive 
understanding of short wave radar reflectivity (Hasselmann et al., 1985). To have a 
more reliable estimate of this parameter, it is crucial to obtain more direct measure- 
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merits of the coherence time which determines the effective azimuthal resolution of 
the SAR and hence INSAR imagery. 

2. Estimation of Radar Scene Coherence Time Based on SAR and 

INSAR Simultaneous Images 

The temporal coherence of the ocean scene can be investigated by using 
simultaneous SAR and INSAR records. As described in section B of this chapter, 
the image intensity of SAR is a map of the covariance of the complex reflectivity. 
The mean image intensity obtained by a single antenna is computed by ensemble 
averaging neighboring pixels, and is directly related to the covariance of complex 
reflectivity: 

( r i(x,<)rJ(x,f + r)> = <7i(x,t)exp [-r 2 /2r c 2 ] (4.52) 

Similarly, for the second antenna separated physically by a distance B from the first 
antenna, the covariance of complex reflectivity is taken to be: 

(r 2 (x, t -f At)rl(x, t -f At + r)) = <r 2 (x, t , At )exp [-t 2 /2t 2 ] (4.53) 

where At = B/2V and a is the radar cross section which represent the covariance 
of short wave reflectivity modulated mainly by long wave tilt. Since the temporal 
variations of these waves are slow compared to the interferometric time interval At, 
it is reasonable to assume that the radar cross-section does not vary during the time 
At. Therefore 

<7i(x,f) ~ c 2 (x,t + At) ~ a(x,f) (4.54) 

When corresponding pixels from the two images are combined to extract 
the interference phase, the temporal cross-correlation between the corresponding 
pixel scattering elements introduces an additional cross-section modulation due to 
the random velocity spread of subresolution intermediate waves modulating the 
short wave reflectivity. This relatively fast time decorrelation of the scene in a given 
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time interval At = B/2V attenuates the INSAR signal (see Equation -(4.33)). The 
magnitude of the cross-product of the complex images is directly related to the 
cross-covariance of complex reflectivity, and is given by: 

( r i(x,t)rj(x,t + Af + t)) = <7j l2 (x,t + At) exp [-r 2 /2r c 2 ] (4.55) 

where 

<ri )2 (x, t, At) = cr(x, l) exp [—At 2 /2r 2 ] (4.56) 

Substituting Equation (4.56) into Equation (4.55) yields: 

(r : (x,t)rj(x,t + At + r)) = a(x,i)exp [—At 2 /2r^j exp [-r 2 /2r 2 j (4.57) 

Neglecting for a moment the additive noise effects, the effective radar 
cross-section obtained by INSAR (Equation 4.56) normalized to the product of the 
square root of both corresponding SAR radar cross-section by invoking Equation 
(4.54) yields 

K _ (n(x,t)r^(x,t + At)) _ 

( r i(x,t)rj(x,t))V2 (r 2 ( x ,f)r 2 (x, t)) 1 / 2 

= exp[—At 2 /2r 2 ] (4.58) 

with 

At = B/2V 

Equation (4.58) expresses the ratio between the INSAR ocean surface cross-covariance 
of complex reflectivity and the two individual SAR covariance of complex reflectiv¬ 
ities. 

Coherent radars such as SAR are inherently subject to speckle or fading 
noise. By averaging neighboring pixels the speckle noise is reduced at the expense 
of nominal resolution. This operation is similar to the multiple look-operation, in 
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which the processor optimal (maximum) integration time is divided into shorter 
integration periods, and by averaging similar images with a time shift, the fading 
noise is reduced at the expense of resolution. 

In the data analysis, relatively large identical areas (> 10 3 pixels) for the 
two SAR and INSAR intensity images were chosen for averaging. The noise power 
is reduced by: 



where N is the number of averaged pixels. The speckle and system noise are further 
reduced by normalizing the ratio K in Equation (4.58) by a similar ratio K 0 obtained 
from a stationary scene, like land, which fortunately was in the images. The scene 
coherency of land varies between fields with objects like plants and trees moving 
with the wind, and urban regions, which should be absolutely stationary and hence 
should ideally have infinite scene coherence time. Due to speckle noise and system 
errors the apparent coherency of urban regions is not infinite; thus normalizing 
K from Equation (4.58) with K 0 obtained from the urban regions of Marina, will 
further reduce the noise resulting in a better estimation of r c using this method. 
Based on this assumption, Equation (4.58) is modified to be: 

K/K 0 = exp[-At 2 /2r 2 ] (4.60) 


and the scene coherence is taken to be 


_ At 

y'-21n(/</K 0 ) 


(4.61) 


Equation (4.61) demonstrates the sensitivity of the scene coherence time estimation. 
A change in the ratio K/Ko of 2% will change the estimation in the scene coherence 
time by almost 20% (see Figure 4.4). The estimates of r c using the above experi¬ 
mental method yields a more credible L-band scene coherence time 0 (100 msec), 
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Figure 4.4: Dependence of scene coherence time on K / Ko- 

which is consistent with the few earlier known estimates based on observations (e.g., 
Plant and Keller, 1983). 
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V. REFRACTION SHOALING AND 
DIRECTIONAL SPECTRA OF NEARSHORE 

OCEAN WAVES 

A. INTRODUCTION 

This chapter presents a brief review of linear ocean surface waves theory in 
shallow waters. A numerical refraction modei (Dobson, 1967), followed by wave 
spectral analysis and measurement methods, are described. The behavior, analysis, 
and in situ observation of ocean surface waves are well developed subjects. How¬ 
ever, many of the concepts presented are difficult to interpret in terms of coherent 
and noncoherent radar ocean surface imagery. The content of the present chap¬ 
ter is essential for interpretation of the experimental results and interdisciplinary 
understanding of the issues involved. 

B. TRANSFORMATION OF NEARSHORE OCEAN 

WAVES 

When waves propagate from deep to shallow water they are transformed as 
they begin to ‘feel’ the bottom. The transformation includes shoaling, refraction, 
diffraction, reflection, and dissipation. In the present analysis only shoaling and 
refraction phenomena are considered, since we ignore reflection (beach slope less 
than 1/20) and assume energy is conserved at steady state and is not translated 
lateraly along the wave crest (ignoring dissipation and difraction). The experimental 
site at Marina provides an opportunity to study wave refraction and shoaling as 
imaged by INSAR and ground truthed by a shallow water array off Marina. 
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1. Wave Refraction and Shoaling 


Starting with the wave energy flux balance assuming an idealized situa¬ 
tion in which wave energy is conserved at steady state and not transferred lateraly, 
yields the simplified energy flux equation: 

F 

(5.1) 


^k + ^ = 0 

dx T dy 


where F (x,y) is the wave energy flux and is constant. Given a monochromatic 
deep water gravity wave field propagating shoreward as depicted in Figure 5.1, the 
following energy flux balance applies 

EoCgobo — E\Cgib\ (5.2) 


where E 0 , E x are wave energy Cg 0 , Cg\ are group velocities and b 0 , b x are distances 
between wave rays in deep and shallow water respectively. We further discurs the 
derivation of a wave height relationship for shoaling waves from Equation (5.2) based 
on linear wave theory. 

2. Refraction 

Assuming the gravity waves are linear and dispersive yields the relation¬ 
ship. 

u> 2 = < 7 fctanh kh (5.3) 


where u> is the wave frequency, g is the acceleration of gravity, k = 27r/L is the 
wavenumber magnitude (L is the wavelength) and h is the water depth. The group 
velocity for a surface gravity wave train is defined as 

Cg ~ dk 

Substituting Equation (5.4) into Equation (5.3) yields the group wave velocity for 
gravity waves 



r _ g (, , 2 “ 

9 2 \ sinh 2 kh 


(5.5) 
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Figure 5.1: Wave refraction and shoaling. 


where C is the wave phase velocity. 

The instantaneous surface elevation of a monochromatic wave is given 
by: 

= asin(k • x — wi) (5.6) 


where a is the wave amplitude and x = (x,y) is the position vector. The variance 
of this wave is 

= jrJ o (^<)-S<)) dx= J ( 5J ) 

where L is the ocean wavelength. Therefore the energy of a monochromatic wave 
in an irrotational and incompressible fluid given by: 

\ P9 a2 = \ PS 1 * 2 ( 5 . 8 ) 


where H = 2 a is the wave height. 
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Substituting Equation (5.8) into Equation (5.2) yields the wave height 
variation equation for shoaling waves (Mei, 1980) 

si 11/2 ri 11/2 

H = Hq = Hq k s kn (5*9) 


where k, is the shoaling coefficient 


(5.10) 


and kn is the refraction coefficient 


**= j • (5.11) 

Within the limits of linear wave theory, the wave frequency u> must re¬ 
main constant as waves propagate from deep water ashore. Using this approximation 
in the dispersion relationship rewritten in a wavelength form results in 


^‘(ir)fir 


(5.12) 


Equation (5.12) illustrates the variation in wave length with depth for shoaling waves 
to preserve the invariant wave frequency. For deep water waves 


(5.13) 


and Equation (5.13) can be rewritten as (Shuchman and Kasischke, 1981) 


L\ — Lq tgh 


(5.14) 


For straight and parallel of bathymetric contours, a simplified refraction 
formula based on Snell’s low can be stated as: 


C 1 = C 0 ^i 

sinao 


(5.15) 
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where C = u>/k is the wave phase speed. Using Equation (5.3) to compute the phase 
speed and substituting the result in Equation (5.15) results in the refraction angle 
a between the wave rays and bottom contours for waves propagating from deep to 


shallow water 


aj = sin -1 sin ao tan h 


(¥)] 


(5.16) 


In the data analysis described in Chapter VII, the shoaling effect is 
estimated by using Equations (5.12) and (5.14), and the refraction effect is estimated 
using Equation (5.16) for straight and parallel contours only. An attempt to estimate 
refraction and shoaling of coastal gravity waves using SEASAT SAR data was done 
by Shuchman and Kasischke (1975). The method presented in this study is seen to 
be more accurate. 

3. The Ray Trace Model 

For complex bathymetric contours which are not straight and parallel to 
the shoreline, a linear numerical refraction and shoaling model by Dobson (1967) is 
used in our analysis. Following the approach established by Munk (1952) on wave 
refraction, Dobson developed a linear refraction model using the characteristics of 
the propagation of curving wave rays. Considering the geometry illustrated in Figure 
5.2 the ray curvature equation in Cartesian coordinates is given by: 


n{x ' s ’ t) = CM [ sma{t) —dT~~ msa{t) ~W~ 

and the ray separation factor ft is stated by Dobson to be: 

^t + pM f+ «(«)/? = o 


(5.17) 


(5.18) 


p(x,y,t) = 2 cos or 


dC(x,y) 


+ sin a 


dC(x,y) 


(5.19) 
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Figure 5.2: Refraction geometry and notation for Dobson’s model. 


, J. 2 d 2 C(x,y) 0 . d 2 C(x,y) , _ 2 d 2 C(*,y)l 

q{x,y,t) = C(x,y) sin a—-2 sin a cos a ^ -bcos a——— 

(5.20) 

The wave ray path defined by R and a is computed with Equations (17)—(20), 
provided the wave phase speed C(x,y) is known. C(x,y) is constructed grid to grid 
over the ocean bathymetry as propagating to the required position. For more details 
of this method see reference (Dobson, 1967). 

In this study, the bathymetry used in Dobson’s refraction model is based 
on data collected originally by NOAA and projected onto a six second (lat, long) 
Modified Universal Transverse Mercator (UTM) grid. The data was screened and 
interpolated using linear triangulation to fill missing points onto a 200 meter rect¬ 


angular grid (Kalid, 1989). 
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C. OBSERVATION METHODS OF OCEAN WAVE 
SPECTRA 

The relation between the distribution of energy in sea surface waves and the 
direction of propagation is of great interest from the forecasting point of view and 
the understanding the generation processes. For the case presented in this study, 
shoaling wave transformations such as diffraction , refraction and reflection (see 
next section) are greatly influenced by the directional spectral characteristics of the 
waves. However, measurements of directional wave spectra requires greater effort 
than that of frequency power spectra. While the latter can be obtained from a wave 
record at a single point, the former requires simultaneous recording of several wave 
components. 

The measuring techniques which are known to be applied can be classified 
into two categories (Goda, 1985) 

1) In situ measurement method 

• Wave gauge array 

• Directional buoy 

•> Two axis current meter 

2) Remote sensing method 

• Optical technique (sterophotogrametry, holography) 

• Microwave technique 

In Chapter VII, frequency and directional spectra obtained from a wave gauge 
array are compared with two-dimensional wavenumber spectra obtained by a remote 
sensing interferometric SAR. Therefore, in the present section a synopsis of the basic 
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theoretical and observational background of direct and remote sensing methods are 
described. 

1. Ocean Surface Wave Field Spectra 

The characteristics of wave field properties can be determined from the 
observed sea surface elevation and the irrotational surface velocity. From linear wave 
theory there is a direct relationship between the sea surface elevation, velocity and 
pressure spectrum. The horizontal component of the wave orbital velocity can be 
described by a product of the ocean surface elevation r/(x, t) and a transfer function 


II „(w): 

tt(0 = Hu{u) r/(x,<) 

(5.21) 

where 

?/(x, i) = a cos(k • x — w<) 

(5.22) 

and 

. cosh k(h + z) 

// ” ( " ) = w smh U ’ 

(5.23) 

where z 

is the depth of observation. The wave amplitude a 

is assumed to be small 

compared with wavelength and depth of water as required by linear theory. Similarly, 

the vertical velocity component is given by 



w(/) = »?(x,0 

(5.24) 

where 

tj / \ . sinh k(h + z) 

//»=-.* sinhM 

(5.25) 

The pressure wave form at depth z is taken to be 



P(t) = Hp(u>) rj(x,t ) 

(5.26) 

where 

„ . , cosh k(h + z) 

cosh kh 

(5.27) 
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In the direct method, time series records are acquired from the random 
ocean processes. Assuming a stationary ocean surface, the classical analysis ap¬ 
proach is the spectral technique, i.e., linear superposition of an infinite number of 
sinusoidal components with small amplitude and slope. Considering the water col¬ 
umn as a constant parameter linear system, the input-output frequency is invariant. 
For input surface elevation spectra, Equations (5.21), (5.24), and (5.26) can be ex¬ 
pressed in terms of input output power spectra related by the corresponding linear 
frequency invariant transfer function H(u): 

S.,,H = |«(*>)PS.(w) (5.28) 

Equation (5.28) reflects the fundamental idea of measuring the ocean wave frequency 
power spectra. The random sea surface elevation is considered as the input and the 
pressure or velocities as output. Next we discuss relevant statistical characteristics 
of the sea surface elevation related to the ocean wave spectra. 

The random ocean surface is characterized by the autocovariance, p( r, r), 
of the surface displacement, 77 as: 

p(r, r) = (t?(x, t) »?(x, r + r, t + r)) (5.29) 

where r is the displacement vector. 

The ensemble average of the autocovariance surface elevation in Equa¬ 
tion (5.29) defines the decorrelation time and space scales of the sea surface. The 
Wiener-Khinchin theorem states that for any ergodic stationary random process 
(e.g., ? 7 (x, t)) with zero mean, the wavenumber-frequency spectrum is given by 
the Fourier transform of the autocovariance function p(r, r) of the random process 
j?(x,t). 

S(k,u>) = —73 J J p( r > T ) exp[-i(k-r-u>r)]dr dr (5.30) 
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with the inverse Fourier transform 


p(r, r) = J J 5'(k,u;) exp t(k • r — wr)jdk du 


(5.31) 


The mean square of surface elevation (ij 2 (t,t)) is obtained from Equation (5.30) for 
r = 0 and r = 0. This variance term is proportional to the wave energy. 

Under the above approximations the Fourier transform of the autoco¬ 
variance of the surface elevation is computed as a function of wavenumber and 
frequency, with position and time as parameters. 

2. Relation to Measurement Methods 

The computation of the frequency wavenumber spectra requires simul¬ 
taneous measurements of surface displacement in space and time. This requires an 
array of wave sensors and an enormous data processing effort. The omni-directional 
frequency spectrum is obtained with less effort and therefore more common, since 
it requires only a single wave sensor at a point, collecting temporal data with no 
spatial information about the wave field. At a given position r = 0, the frequency 
spectra is given by 

S(u) = / 5(k,w) <fk = J- //>(0, r) exp [iu>r]dT (5.32) 


Since only positive circular frequencies have physical meaning and S(io) is symmetric 
about u = 0, Equation (5.32) can be rewritten as: 

1 r+°° 

S(u>) = — / p(r) cost or dr (5.33) 

tt Jo 

Since the total energy for a given sea state is unique, the integration areas 
of the frequency and wavenumber spectra are related by the Jacobian, yielding the 
unidirectional wavenumber spectrum 

S(k) = ^ S{ u.) (5.34) 
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where k = |k| and du>/dk = C g the group wave velocity. Using the linear dispersion 
relation u = u(k) in computing C g ,S(k) is obtained. This transformation can not 
be used in the high wavenumber, short waves part of the spectrum, since the intrinsic 
frequency of the shortest waves w 0 (fc) are Doppler shifted due to modulation by the 
underlying dominant waves and/or currents on which the short waves are borne. 
Therefore the circular frequency of the short waves is given by 

u(k) - u) 0 (k) -f k • (U 0 b(z> t) + UcnOr, t)) (5.35) 

wheie Uob(x,2) is the slowly varying orbital velocity and U cn{x,t) could be any 
underlying varying currents. Equation (5.35) illustrates the nonuniqueness of the 
relationship between frequency and wavenumber. 

For high resolution microwave radar sensors the backscattering energy 
is determined by the short scale (cm) surface roughness which is borne on an ad- 
vecting surface. Equation (5.35) replaces the linear dispersion relationship between 
frequency and wavenumber (Equation (5.3)). 

3. Wave Directional Spectrum 

As noted in the literature review (Equation (2.7)) the normalized radar 
cross-section cr 0 >s directly proportional to the two dimensional wavenumber spectrum. 

V>(k) = J S{k,u)<ko (5.36) 

The instantaneous autocovariance can be expressed as: 

p(r,0) = J J ip(k exp[-i(k • r)]dk (5.37) 

and inversely 

V»( k ) = ^j~)2 / ^ ex P[ { ( k * r )\ dr ( 5 - 38 ) 
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To obtain V’(k), one needs spatial observations of the surface elevation autocovari¬ 
ance (Equation (5.38)). Attempts to measure this phenomena are still under exten¬ 
sive research using remote sensing techniques such as sterophotogrametry (Shemdin, 
et al., 1988) laser optics (Keller and Gotwols, 1983) or Doppler radars (Plant and 
Keller, 1983). 

The two-dimensional wavenumber spectrum is also claimed to be com¬ 
puted from SAR images. However, there is still doubt about the precise mechanisms 
whereby waves are imaged by SAR. Therefore, the question remains whether the 
SAR image intensity spectrum bears much resemblance to the directional wave spec¬ 
trum. In this study, a new approach is suggested for obtaining a two-dimensional 
wavenumber spectrum using interferometric SAR imagery. This method has a 
clearly observed imaging mechanism based on the measurement of the varying or¬ 
bital velocity Uob{x,L) of the dominant surface ocean waves. The two-dimensional 
wave number spectrum of INSAR’s phase image represents the spatial distribution 
of the surface velocity power in the viewing direction. In describing this spectrum, 
assume for the moment a monochromatic dominant wave oriented at an arbitrary 
angle (f> relative to the flight direction (see Figure 4.1) and an incidence viewing 
angle 7 . The surface elevation of this wave as seen by the imaging platform is given 
by: 

t/(x, t) = a g{~f,(f>) cos(k*x — c ot) (5.39) 

where a is the wave amplitude and 

< 7 ( 7 , (f >) = (cos 2 7 + sin 2 7 sin 2 <^) 1 ^ 2 (5.40) 

After removing the current and phase velocities, the residual imaged 
radial velocity component UoBr can be computed. Using deep water linear wave 
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theory 


UqBt (x) = — gp = aw g(r, <fi) sin(k -x-wt) 


(5.41) 


Then the two-dimensional wavenumber spectra for the orbital velocity is 

1 r+L/2 r+L/2 

^Uobtv^zi ^v) = I I ^OBri^^y') 6xp[ i[k x x -f- fcyt/)] dxdy (5.42) 
L‘ J-L/2 J-L/2 

The spectrum described in Equation (5.42) is tested against directional spectra 
obtained from a wave pressure array. A synopsis review on the theory of in situ 
measurement of directional spectra is described next. 

Longuet-Higgins, Cartwright and Smith (1963) developed a directional 
wave bouy based on measuring the angles of pitch and roll and the heaving accel¬ 
eration of the ocean surface. The heave is directly related to the surface elevation, 
while the pitch and roll are related to the two horizontal orthogonal components of 
wave surface slope namely, dr\/ox, d'j/dy. The directional spectrum S(f,0) can be 
expressed as a product of the one-dimensional frequency power spectrum S(f) and 
an angular spreading function £)(/, a) 

S(f,a) = S(f) D(f,a) (5.43) 

where a is the wave direction. D(f,a ) is associated with the distribution of energy 
in the frequency power spectrum, over direction and has a unit area over direction, 

f D(/,a) da = 1 (5.44) 

Jo 

An objective of measurements is to determine the characteristics and criteria in 
selecting £>(/, a) which most closely match the observed angular distribution of 
wave energy density. Using the most common in situ measurement technique - a 
floating buoy or slope array of four sensors in a square, one has to address their 
limited spatial capabilities in defining directionality, namely £>(/,a). 
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Following Longuet-IIiggins et al. (1963), the directional spreading func¬ 
tion £>(/, a) can be expanded as a Fourier series: 

D(f,a) = i Q + £[a n cos(na„) + b n sin(na„)]|j (5.45) 

For the previously discussed bouys or slope arrays which measure only the surface 
slope and displacement, only the first two harmonics can be estimated. The hori¬ 
zontal sea surface slopes are obtained by taking the corresponding spatial derivative 
of the sea surface elevation (Equation 5.6) 

*?r(0 = = (—t&cos a) T)(t) = H x (f ) t ]{t) (5.46) 

Vy(t) = = (“** cos «) = Hy if) V(0 (5.47) 

Five Fourier coefficients can be determined from the co(C) and quad ( Q ) spectra of 
the surface elevation and surface slope. The spectrum of r/ r (f) is given as 

C„„U)= r\H,(!)\'S(S,«)da (5.48) 

Jo 

Substituting Equations (5.43), (5.46) into Equation (5.48) yields 

SwM) ~ f k 2 cos 2 aS{f)D(f,a) da = k 2 S{f) f cos 2 aD(f,a) da (5.49) 
Jo Jo 

Similarly the rest of cross-spectral quantities are obtained 


C„,(f) = k 2 S(f) r sm* a D(f, a) da 

Jo 

(5.50) 

C„(J) = S(f) pD(/,a)da=S(f) 

Jo 

(5.51) 

y(/) = k 2 S(f) 1 cos a sin aD(f,a) da 

JO 

(5.52) 

<?,„(/) = kS(f)J\osa D(!,a) da 

(5.53) 

Qnvjf) ~ kS (f) [ sina D (fi a ) da 

J 0 

(5.54) 
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From Equations (5.49)-(5.54) Longuel-Higgins et al. (1963) computed the Fourier 


coefficients of £>(/, a) as follows: 

a 0 = — [ D(f,a) da — - (5.55) 

'K JO 7T 

ai = -[ D(f,a) cos a da = QvnM) (5.56) 

7T JO t TCK 

bi = ~J 0 W> a ) sin a da - Qw v (f) ( 5 - 57 ) 

\ J o 2n D(f, a) cos2a da = ^ (C, tt , x {f) ~ C„ yVy (f)) (5.58) 

bi = ~J Q D (f' tt ) S ‘ n 2a da = ^2 ( 5>59 ) 


Longuet-Higgins et al. (1963) used these five Fourier coefficients to de¬ 
scribe an assumed cosine bell angular energy distribution function. This method 
can only describe a unimodal energy distribution with limited resolution, which is 
frequently not realistic. These coefficients can be represented by an infinite number 
of angular spreading functions. It is therefore desirable to seek a technique (if pos¬ 
sible) to match the coefficients obtained from the data to a more realistic angular 
distribution function. 

Grauzinis (1990) improved the Longuet-Higgins et al. (1963) method 
in matching analytically a more realistic angular distribution function D(a) to the 
obtained Fourier coefficients. The technique is based on representing the Fourier 
coefficients of D(a) as the superposition of a weighted sum of symmetric unit area 
density functions oriented at different direction 

a it + ibk = X)wj m k {j) exp {jkfy) (5.60) 

j 

where w ; - is the weighting parameter, m k is the canonical coefficient of the even unity 
spreading function and (3j is the angular orientation of component j. For a typical 
pressure array (four sensors) or pitch and roll buoy, the obtained cross spectral wave 
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quantities result in the Fourier coefficients of Equation (5.60) through the second 
order and are taken to be 


Ji + ibi = mi [w, exp (ifa) -f w 2 exp(t&)] 

(5.61) 

a 2 -f tf> 2 = m 2 [wi exp(*2/?i) + w 2 exp(z‘2/? 2 )j 

(5.62) 


To solve the system represented by Equations (5.61), (5.62), further as¬ 
sumptions have to be made. The most sufficient assumption made by Grauzinis 
(1990) is that the two canonical coefficients are related by a power law r 

m 2 = (5.63) 

where r defines the broadening (r increase) or the sharpness (r decrease) of the 
distribution function defined initially by the canonical coefficients. For a unimodal 
wave field and m = 0.9 and r = 4, the half power distribution is 12 degrees, which 
is about one third of the Longuet-Higgins et al. (1963) cosine bell resolution. The 
improvement of this technique relative to the Longuet-Higgins et al. (1963) method 
is the variability in the sharpness of the model distribution function, which enables 
it to be matched with better directionality to the observed energy distribution func¬ 
tion. Also. Equations (5.61), (5.62) introduce two distributions of wave energy with 
different orientations and weights. This allows the model to represent a bimodal 
wave system which is not possible by the Longuet-Higgins et al. (1963) method. 

This exact method of matching the Fourier coefficients to a wave spread¬ 
ing function is not unique. The question that remains is how to optimize the selec¬ 
tion of the distribution that best matches the real data. Cartwright (1963) suggests 
that a criteria to improve the selection of the distribution function can be obtained 
from the ratio of the magnitude of the Fourier coeffients. Grauzinis (1990) further 
suggests a beam separation index criteria to indicate the best distribution function 
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fitting. This criteria aids the selection of the canonical coefficients which fulfills the 
analytic solution to give zero beam separation from the required distribution. For 
different r’s, the bimodal distribution with the minimum peak separation is selected 
as the best fit. Despite the analytical efforts to improve the directional resolving 
wave power obtained from pitch and roll buoys or pressure/current limited arrays, 
an rms resolution of 4-6 degrees appears to be the best that can be achieved using 
this method. 
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VI. THE MONTEREY BAY WAVE FIELD 
AND CURRENTS ‘89 REMOTE SENSING 

EXPERIMENT 

A. INTRODUCTION 

The primary objective of the experimental part of this study is to verify the 
feasibility of an airborne INSAR method to reliably spatially sense ocean surface 
wave and current fields. The verification is done by comparing the results of this 
unique technique with simultaneous in situ conventional temporal and local spatial 
measurements including directional wave spectra from an array, local wind and sur¬ 
face drifters measurements. To simplify the interpretations, the time of the year 
and the site of the experiment were selected so as to provide near linear condi¬ 
tions of wave imaging (narrow band wave spectrum of low to moderate variance) 
and a moderate sloping beach (approximately 20 m/km) allowing sufficient locally 
homogeneous conditions for linear spectral analysis. 

The relatively shallow water provides an opportunity to study for the first 
time wave refraction and shoaling using INSAR. Four partially overlaped areas (~ 
12 km x 6 km) were imaged by INSAR from four different directions centered on 
the shallow water wave array off Marina. Results from three legs are reported in 
this study. 

Overall, three field tests were performed in the vicinity of Monterey Bay: 

• On June 1, 1989, several prototypes of passive radar drifters were 'ted at 
sea to measure ocean surface currents using Lagrangian bouys. The drifters 
were designed and built at the Naval PostGraduate School and were tested 
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from the Research Vessel (R/V) Ricketts in the vicinity of Moss Landing (see 
Appendix C). 

• On July 23-24, 1989, a preliminary coordinated simultaneous airborne and 
ground truthing experiment was performed. The experimental scene was cen¬ 
tered on the wave and wind stress NDBC pitch and roll buoy off Monterey Bay 
(see Figure 6.2). NASA/JPL DC-8 radar overflights were conducted on July 
24,1989 between 1000 and 1040 Pacific Day Time (PDT). The R/V Point Sur 
on a cruise from 232200 to 242400 July 1989 acquired ground truthing data 
for comparison with the following near ocean surface phenomena: 

- Vertical profile of current vectors within the mixed layer measured by 
Acoustic Doppler Current Profile (ADCP). 

- Surface currents sensed by Lagrangian floating buoys which were posi¬ 
tioned by Loran-C and Mini-Ranger navigation systems (see Section E4). 

- Wind speed and direction which generates the short ocean surface wave 
field patterns responsible for the resonant radar backscatter (Bragg waves). 

— Vertical density gradients at the depth between 10 and 250 m measured 
by Conductivity Temperature Depth (CTD), which might affect surface 
dynamical processes. 

Wind stress and dominant wave time series records were acquired from the NOAA 
pitch and roll buoy. Unfortunately, the Interferromatric SAR images could not be 
obtained from this experiment. An error in setting up the fore and aft radar anten¬ 
nas with different polarizations (vertical and horizontal), instead of being identical 
(both vertical or horizontal), prevented fringes formation essential for creating an 
interferogram. 


76 







Based on these preliminary experiments, the experimental plan was modified 
for a following experiment. The field test site was moved to an area close to the shore 
to include stationary land in the INSAR imaging scene to assist aircraft attitude 
corrections in addition to the aircraft Inertial Navigation System (INS) (see Chapter 
III). A second modification was to replace the passive radar drifters with thin floating 
plywood sheets (see section F3) to minimize windage effects. The following and final 
experiment for this study was carried out on September 8, 1989 in Monterey Bay 
centered on the shallow water wave array off Marina Beach. An illustration of the 
July and September experimental scenarios is depicted in Figure 6.1. A detailed 
description of the Marina experimental site is given in section B, followed by a 
description of the Marina experiment in section C. The environmental conditions 
for the experiment are given in section D. 

B. THE EXPERIMENTAL SITE 

The Monterey Bay is characterized by variable bathymetric features including 
the largest submarine canyon in the western hemisphere bisecting the Bay (Figure 
6.2). These bottom features cause substantial perturbations on the near shore wave 
fields and tidal currents. The experimental site for data acquired during this study 
is in the nearshore region of Marina within Monterey Bay. This region is shadowed 
by Point Santa Cruz in the northwest and Point Pinos in the southwest. 

The Marina site was selected because of the location of a wave array for 
ground truthing and the required inclusions of land in the scene. Also, the site 
has a relatively continuous, though moderate, alongshore bottom slope (about 20 
m/km) resulting in a relatively small but sufficient region of locally homogeneous 
bathymetry. This fact constrained the spectral analysis to be applied to small 
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spatial regions to meet the required homogeneous conditions. This area presents an 
interesting problem in predicting directional wave spectra. 

The sloping bottom introduces additional various dynamic ocean features 
like gradients in the tidal flows. The present study concentrates on the dominant 
surface waves spectra and currents. Shoaling waves undergo refraction, reflection 
and diffraction influenced by the wave approach angle and bathymetry. The complex 
bathymetry of Monterey Bay is shown in Figure 6.2, and the detailed bathymetry 
near Marina in Figure 6.3. Wave refraction at the Marina array site is shown in 
Figure 6.4 for different frequencies (0.03 Hz to 0.13 Hz). Linear refraction theory 
(see Chapter V) predicts that lower frequency waves (longer waves) experience more 
refraction. All possible incoming ray traces with wave angles at increments of 0.2° 
using Dobson’s model (see Dobson, 1967) are plotted from shallow to deep water. 
High frequency waves of at least 0.11 Hz are limited to directions from 234°-311°, 
whereas lower frequency waves may arrive aiso from southerly directions since they 
experience more refraction. 

An important feature in the experimental site is the Salinas River, which was 
blocked from the sea at this time of the year (June-September) providing a natural 
closed basin. This allowed the study of Bragg backscattering waves only, without 
the contamination by interactions with long waves or any underlying currents. 

C. THE MARINA EXPERIMENT 

The Marina experiment consisted of four SAR and INSAR overflights in the 
interferometric mode on September 8, 1989. Radar data were acquired about an 
area centered on the shallow water wave array (see Figure 6.5) with the antenna 
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Figure 6.3: Detailed bathymetry 
Beach. 
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Figure 6.4: Refraction wave ray traces 
frequencies. 
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Figure 6.5: Flight pattern and scene locations for radar data collection 
on September 8, 1989 off Marina Beach. 
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beams aimed at the in situ wave array with a depression angle of 45° and mean alti¬ 
tude of 8250 m. The radar overflights were conducted around 1300 Pacific Daylight 
Savings Time close to the time of maximum tidal flood current. The flight pattern 
and the scene locations of radar data collection are illustrated in Figure 6.5. Some 
key system and flight characteristics are depicted in Table 6.1. More detailed SAR 
system characteristics are given in the next chapter. 

The Lagrangian floating drifters for tracking the very top surface currents 
were deployed and tracked by the R/V Ricketts with the aid of a Piper Cherokee 
Arrow (PCA) light aircraft (Figure 4.15). A portable LORAN-C enabled the light 
aircraft to locate roughly the drifters while the vessel, and hence the drifter positions, 
were recorded using a high accuracy Motorola Mini-Ranger Falcon 484 system. 

The actual time schedule and sequence of events in the Marina experiment 
on September 8, 1989 were as follows: 

081000 Departure of R/V Ricketts from Moss Landing home base. 

081050 On station (Marina shallow water array) to test the navigation system 
(Mini-Ranger) and to estimate local surface currents by deploying and 
positioning two drifters. This pre-measurement was done in order to op¬ 
timize the deployment pattern of the drifters to obtain maximum coverage 
of variable surface currents within the expected airborne radar scene for 
data collection. The variability of the bathymetry was part of the concern 
in ground truthing the radar surface current estimates at variable depth. 
081120 Termination of initial test, vessel steaming westward to the first deep 
water deployment station. 

081145 - 081230 Buoy deployment from west to east in three distinct groups, 
each group consisting of five drifters for a total of fifteen (see Figure 
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TABLE 6.1: System and Flight Data for the Different Aerial DC-8 Air¬ 
craft Coverage During Marina Experiment, September 8, 1989 



90° 

EAST¬ 

WARD 

LEG 

270° 

WEST¬ 

WARD 

LEG 

360° 

NORTH¬ 

WARD 

LEG 

Initial data collection 
time, PDT 

12:40 

13:01 

12:51 

Aircraft velocity, m/s 

214 

216 

214 

Aircraft altitude, m 

8285 

8273 

8178 

Pulse repetition frequency , Hz 

283 

285 

282 

Near incidence angle, degrees 

17 

Nadir 

Nadir 

Far incidence angle, degrees 

52.7 

50.2 

50.7 

Ground pixel range spacing, m 

8 

7 

6 

Azimuth pixel spacing, m 

12.1 

12.1 

12.1 

Effective coverage 
range x azimuth, km 2 

6 x 14 

5.1 x 14 

5.8 x 14 


6.6). Group one drifters, numbered 1 thru 5, were deployed furthest 
ofTshore approximately 3.4 nm west of the shore where the mean water 
depth is 71 m. This group was positioned four times during the three 
hour observation period. Group two drifters, 6 thru 10, were deployed 
closer to the shore, approximately 1.6 nm west of the beach where the 
mean water depth is 44 m. These were positioned three times. Group 
three drifters, 11 thru 15, were deployed less than 0.6 nm from the beach 
and washed ashore before a second position could be obtained. Within 
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each group, the drifters were deployed roughly 0.1 nm apart as the vessel 
steamed eastward toward the beach. 

081220 The Piper Cherokee Arrow light aircraft arrived on station Monterey Air¬ 
port and started to position and photograph the drifters from west to 
east. 

081230 - 081515 Drifters chased, with vessel heading west to position the drifters 
first deployed. The initial detection of the Lagrangian drifters was quite 
tedious. Applying Direction and Range (DR) estimates based on the 
early test results, the drifters were not detected. The light aircraft was 
very helpful in detecting the drifters initially, although the radio com¬ 
munications broke down. Therefore the actual communication was done 
by turning on the landing lights and heading towards the R/V Ricketts 
once the Piper overflew a group of drifters. As predicted, the onshore 
tidal current became more intense close to the time of midtide range 
(1320), yielding a stronger eastward component, which was measured by 
the drifters. The second DR prediction based on the initial physical po¬ 
sitioning was close to the actual positions of the drifters. 

081240 - 081310 SAR overflight acquired data in the INSAR mode. 

081405 The Piper light aircraft heads home to Monterey. 

081515 R/V Ricketts heads home to Moss Landing after retrieving ten of the 
drifters from the two most western groups. The most eastward group of 
drifters, which were washed ashore, were collected four days later from 
Marina Beach. 
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Figure 6.6: Drifter deployment positions and Mini-Ranger reference sta 
tions during the Marina experiment, September 8, 1989. 
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D. ENVIRONMENTAL CONDITIONS 


1. Introduction 

The environmental conditions which occurred during the Marina exper¬ 
iment on September 8,1989 were typical late summer conditions for Monterey Bay. 
Light westerly winds (less than 4 knots), surface air temperature of 14-16°C and 
clear skies up to 1000 feet blocked by a stratus cloud deck, were observed during 
the 1000-1500 Pacific Daylight Time (PDT) experiment period (see Table 6.2). 

The mild local meteorological conditiv ns during this period resulted in 
almost no local generation of sea. The sea surface was predominantly long crested, 
narrow-band waves propagating shoreward. Therefore, the scene was linear in char¬ 
acter, yielding a linear Modulation Transfer Function (MTF) of the ocean surface 
for imaging radars. 

Due to the lack of predicted tidal current information for Marina Beach, 
the time of maximum tidal current was estimated using predicted information of 
tides at Monterey Harbor and Moss Landing, which are located almost in equal 
distance (17 km) south and north off Marina Beach. Predicted tidal high and low 
water times (NOAA 1978 Tide Tables) of Monterey Harbor and Moss Landing are 
plotted in Figure 6.7. Based on the figure, the maximum flood tidal current was 
probably at 1320 PDT. 

The summer climatological atmosphere surface circulation in the eastern 
North Pacific consists of a high east-west pressure ridge stretching from west coast 
California to almost 140° E (east longitude). The ridge extends to approximately 20° 
N (north latitude). Below 20° N, the atmosphere has tropical characteristics. This 
late August high pressure system was deformed in shape and orientation (see Fig¬ 
ure 6.8) by a series of northwest lows, which might be the source of the northwesterly 
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TABLE 6.2: Environmental Conditions During the WAFICIN Marina 
Experiment 


Time and date 

Air temperature 
Surface Water temperature 
Relative humidity 
Wind direction 

Wind speed 
Visibility 

Significant wave height 
Dominant swell periods 

Expected time of maximum 
flooding tidal current 

Swell directions 


1030-1500 PDT, 

September 8, 1989 

14 to 16°C 

~ 13°C (est.) 

~ 60% to 70% (est.) 

Varying between northeast 
and southeast 

< 2 m/s 

Up to 1000 ft 7 miles. 
Above 1000 ft overcast by 
low stratus cloud deck. 

< 0.6 m 

a. 9.1 sec 

b. 15.9 sec 
1320 PDT 

a. East Northeast 
(long swell) 

b. East Southeast 
(short swell) 
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Figure 6.7: Monterey Harbor and Moss Landing predicted tides, Septem¬ 
ber 7-9, 1989. 
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9.1 sec. swell with travel time of almost 4 days observed in the Marina experiment 
on September 8, 1989. The 15.9 sec. long crested westerly swell is postulated to be 
generated by a late winter southern hemisphere storm located at 40° S to 60° S in 
the South Pacific. The travel time for the swell from this storm to Monterey Bay is 
almost 8 days (see Figures 6.8 and 6.10). 

2. Evolution of Environmental Conditions 

When investigating the possible sources of the observed swell waves in 
the Marina experiment, three different sources of process data were analyzed: 

• Surface winds and pressure field patterns between August 31 and September 
8,1989, based on NORAPS/SAC model outputs processed by Fleet Numerical 
Ocean Center (FNOC) in Monterey (see Figures 6.8, 6.9). 

• Wave group velocity vectors and significant wave height contours, based on 
Global Spectral Ocean Wave Model (GSOWN) outputs, between August 31 
and September 8, 1989, also processed by FNOC in Monterey (see Figures 
6 . 10 , 6 . 11 ). 

• Geostationary Operational Environmental Satellite (GOES) visible imagary 
(Figure 6.12). 

3. Surface Winds and Pressure Fields 

The surface pressure field pattern on August 31, 1989 in the eastern 
North Pacific was in good agreement with climatology. From the 1st of September 
1989, the east-west orientation of the ridge axis shifted to a northeast, southwest 
orientation. This shift was caused by the development of a strong low pressure 
system south of the Aleutian Islands and the development of an induced trough 
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Figure 6.8: Mean sea level pressure contours and wind vectors at 1700 
PDT on August 31, 1989. 







80N 


40N 


40S 


80S 


Figure 6.9: Mean sea level pressure contours and wind vectors at 1700 
PDT on September 8, 1989. 









Figure 6.10: GSOWN wave group velocity vectors and significant wave 
height contours at 1700 PDT on August 31, 1989. 
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Figure 6.11: GSOWN wave group velocity vectors and significant wave 
height contours at 1700 PDT on September 8, 1989. 
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over California and the Pacific Northwest. These pressure field conditions were 
maintained through September 8, 1989 yielding stable conditions with light winds 
parallel to the coast and westerly light winds between 35° N to 38° N (Figure 6.9). 
These computed wind fields based on NORAPS/SAC model arc in agreement with 
the observed wind tendancy during the Marina experiment. 

The migration of the pressure ridge northward inducing stable light wind 
conditions over the west coast enabled wave systems produced by northwest Pacific 
storms to propagate eastward without local air-sea disturbances. As mentioned 
previously, these developing storms appear to be the source of the shorter period 
swell (9.1 sec.) observed on September 8 in the Marina experiment. 

4. Satellite Data 

Analysis of GOES imagery confirms the presence of an induced trough 
over California and Pacific Northwest (see Figure 4.12). The low pressure system 
south of the Aleutian Islands, which is believed to be the source of westerly winds, 
caused the 1000 foot ocean stratus deck to approach the central and North California 
coast. The winds were not strong enough to advect the clouds over the land. The 
cloud deck prevented the collection of satellite sea surface temperature data for the 
region of interest on the experiment day, September 8, 1989. 

E. INSTRUMENTATION CHARACTERISTICS 
AND PERFORMANCES 

The relevant performance of assets which aided the collection of data in the 
simultaneous experiments are described in the first part of this section. The proce¬ 
dures of data acquisition and initial sorting are described in the second part of this 
chapter. 
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The sensors which were used in the Marina experiment on September 8,1989 
include remote and in situ assets measuring synoptic spatial variations and Eulerian 
and Lagrangian temporal variations of the ocean surface, respectively. 

1. Platforms 

• NASA/JPL DC-8 airborne INSAR and SAR imaging radar remotely sensing 
ocean surface velocity fields and surface intensity reflectivity. 

• R/V Ricketts for deploying and tracking Lagrangian floating drifters following 
the very near surface currents. High accuracy Mini-Ranger navigation system 
(~ 5m RMS error) for positioning the vessel, and hence drifter positions, 
was installed on board of R/V Ricketts in addition to the inherent Loran-C 
navigation system with lower accuracy (~ 65 m RMS error). 

• Light aircraft Piper Cherokee Arrow (Figure 6.15) for detecting and observ¬ 
ing the wooden sheet floating drifters, which were undetectable initially by 
conventional means (radar, visible) from the surface platform (R/V Ricketts). 

2. In situ Assets 

• Ocean wave directional spectra and surface wind data acquired from: 

- Shallow water pressure gauge array located at a depth of 16 m off Marina 
Beach at: 

36°42' North Latitude 
121°48.9' West Longitude 

- Anemometer on Marina Beach located at 36°42' N. Lat, 121°48.9 A W. 
Long from which directivity of the Bragg waves may be inferred. 

• Lagrangian drifters following the very near surface currents tracked by the 
Motorola Mini-Ranger navigation system. 
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• NAVSTAR Global Positioning System (GPS) extremely accurate positioning 
system for ground truthing the geographic coordinates of SAR and INSAR 
images (see Chapter VII. B3). 


3. Platform Characteristics 
a. NASA/JPL DC-8 

The NASA/JPL DC-8 (Figure 6.13) has a multifrequency, multipo¬ 
larization airborne SAR system operating at L, C and P-Band (Held et al., 1988) 
with a full digital recording capability and an inertial navigator system for measur¬ 
ing three dimensionally the aircraft attitude. Detailed performance of the system 
characteristics are given in in Table 6.3. 

At L and C-Band the radar can operate in an along track interfer¬ 
ometer mode by using the physically separated fore and aft antennas mounted on 
the left looking side of the aircraft fuselage. In this experiment, the L-Band inter¬ 
ferometer was used in vertical polarization since the ocean reflectivity at moderate 
incidence angles 20°-70° (Bragg resonant backscatter) is more intense as compared 
with horizontal polarization. 

The receiver has no Sensitivity Time Control (STC) or Automatic 
Gain Control (AGC) capabilities. Instead a wide dynamic gain range (30-56 dB) 
is implemented to compensate for signal saturation. This is achieved using 8 bit 
Analog to Digital Converters (ADC) which are clocked together to ensure registered 
samples. Due to the high digitizing rate of the ADC’s (45 MHz), real data samples 
(the imaginary part is not recorded) are produced at a total data rate of 20 to 60 
MBytes/second. The data from each receiving channel is multiplexed together to 
form an extended “range line” for recording onto any of the three available High 
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Density Digital Recorders (HDDR) at a rate of 80 MBits/s (see Fig. 6.15). No high 
resolution processing capabilities exist onboard the aircraft (Held et al., 1988). 

b. R/V Ricketts 

R/V Ricketts (see Figure 6.16) is a 35 foot light craft suited for day 
cruises on Monterey Bay. The vessel is owned and operated by Moss Landing Marine 
Laboratories. For her size the R/V Ricketts is reasonably equipped with naviga¬ 
tion and deck handling equipment to function as a field station for oceanographic 
observations. 

4. In Situ Assets Characteristics 

a. Shallow Water Wave Array 

The Marina shallow water wave array consists of four pressure gauges 
configured in a square, six meters on a side, aligned 6° off true North (see Figure 
6.16). The pressure sensors are located 0.5-1 m above a mean bottom depth of 
16.2 m. The northern sensor is a parascientific model 2100AS, digiquartz pres¬ 
sure sensor which produces output voltage signal proportional to the absolute water 
surface elevation. The other three pressure sensors are based on a piezo-resistive 
gauge transducer manufactured by Kulite Semiconductor. The array was developed 
to provide wave directivity in shallow water. The azimuthal resolution of the ar¬ 
ray depends on processing methods, and was 4-6 degrees for the range of resolved 
wave periods between 5-20 sec. Data from the shore station is acquired and ini¬ 
tially screened by a computer based control station located at Scripps Institution of 
Oceanography, La Jolla, California. 

b. Wind Speed and Direction from Skyvane Wind Sensor 

A 2102 Skyvane four bladed, low threshold propeller wind sensor is 
mounted on the beach at Marina at an altitude of 27 m from MSL (see Figure 6.17). 
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Figure 6.13: NASA/JPI/ DC-8 research aircraft. Top: overall view: Bot¬ 
tom: SAR antennas. 
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TABLE 6.3: NASA/JPL DC-8 Airborne Imaging Radar Parameters 
Employed During Marina Experiment, September 8, 1989 


Data 

L-Band 

Frequency, MHz 

1237-1260 

Center frequency, MHz 

1248.75 

Wavelength, cm 

24 

Chirp length, (is 

11.25 

PRF/V, Hz / m/s 

1.32 

Band width, MHz 

20 

Peak power, kW 

6 

Azimuthal beamwidth, deg 

8 

Vertical beamwidth, deg 

44 

Antenna beamcenter gain, dB 

18.3 

Azimuthal pixel spacing, m 

3.03 x 4 

Slant range spacing, m 

6.6 

Polarization 

VV 

Receiver dynamic range (dB) 

30-56 

Number of bits per sample 

8 

Number of looks 

1 

Available number of pixels 
(range x aximuth) 

750 x 1024 
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Figure 6.14: _ NASA/JPL airborne radar and recording block diagram. 
The L-band interferometric mode was used for the acquisition of data in 
this study (from Held et ah, 1988). 



Figure 6.15: Light aircraft Piper Cherokee Arrow participated in the 
Marina experiment, September 8, 1989. 
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Figure 6.16: R/V Ricketts participated on the final experiment on 
September 8, 1989. 

The accuracy of wind speed measurement is ± 0.5 m/sec for wind speeds less than 
15 m/sec. The accuracy for wind direction is ± 1% of the full scale (±180°). 

c. Lagrangian Drifters 

Two types of drifters were tested in our experiments to measure near 
surface currents. In the preliminary experiment on July 24, 1989, a passive radar 
free floating drogue was used. It was tracked by a surface navigation radar installed 
on R/V Point Sur and positioned accurately by the Mini-Ranger navigation system. 
A prototype test of various drogue types were tested on June 1, 1989. For further 
details see Appendix C. 

Based on the results of the preliminary experiments, a Lagrangian 
drifter, made of a 2 cm thick plywood sheet and 0.4 x 1.2 m in dimension was used 
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for ground truthing in the analysis of this study. The dimensions of the drifters 
were selected after considering the measurements by Stewart and Joy (1974) and 
Goldstein et al. (1989), in which radar measurements of near surface currents were 
compared with simultaneous measurement by drifters. At L-band the corresponding 
ocean resonant Bragg waves (21-25 cm at 30° viewing angle) mostly feel currents 
in the upper 2 cm of the ocean surface. Therefore, the drifter’s thickness insured 
measurement of the near surface flows experienced by the Bragg backscattering 
waves while the horizontal dimension reduced the sensitivity to vertical motion 
of short waves (Goldstein ct al., 1989). The drifters were painted red with an 
identification number on both sides in case of rolling over in the surf zone (see 
Figure 6.19). 

d. Mini-Ranger Navigation System 

Vessel positions, and hence drifter positions, were determined using 
the Motorola Mini-Ranger Falcon 484 system. The system is widely used for precise 
positioning in a variety of nearshore marine applications. The system consists of a 
transceiving master station borne on a vessel, and up to four stationary reference 
stations (see Table 6.4) along the shore. The range to each reference station is com¬ 
puted at the master station based on the measured round trip travel time of the 
microwave transmitted pulse. The system operates on the principle of a pulse coded 
transponder. The master and reference stations transmit on different frequencies in 
the microwave band, 5510 and 5480 MHz to avoid interference. The standard system 
has a maximum range of 37 km (20 nm) with reference station gain antenna of 13 
dB, and a master station omni-directional antenna gain of 6 dB. Ranges from up to 
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Figure 6.17: Marina shallow water wave array. Top : hardware configu¬ 
ration; Bottom : geometric orientation. 
















Figure 6.19 Plywood drifter used to Lagrangianly follow surface cur¬ 
rents during the Marina experiment, September 8, 1989. Top : perspec¬ 
tive view; Bottom', apparent drifter. 
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TABLE 6.4: Station and Equipment Data 


Horizontal Datum: North American Datum of 1983 (NAD 83) 


Computational Grid: Universal Transverse Mercator (UTM) 

Zone 10, central meridian 123° W 

Grid Scale Factor 

Used in Computations: 0.999743 


Source of Station Surveys by NPS Hydrographic Sciences 

Positions: • 1988-1989 


Reference Station Positions 
(see Figure 6.6): 


Nick- 

Name 

N Latitude 

W Longitude 
(deg min sec) 

UTM Coordinates 
Easting (m) 
Northing (m) 

Elevation 
MSL (m) 

WATS 

36 54 26.895 

602947.7 

23.1 


121 50 39.819 

4085231.4 


HAYS 

36 38 33.826 

607621.3 

137.2 


121 47 45.895 

4055915.3 


DOP 

36 36 5.601 

600434.0 

10.0 

RM1 

121 52 37.479 

4051260.0 


MBAR1 

36 37 6.868 

597861.5 

30.7 

RM 

121 54 20.156 

4053118.5 
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four reference stations, with an rrns triangulation position of ~ 5 m, were recorded 
onboard (see next section). 

F. DATA ACQUISITION AND INITIAL DIGITAL 
PROCESSING 

In this section, the initial processing methods of the data acquired in the 
Marina experiment on September 8, 1989 are discussed. The processing of the 
airborne radar data at JPL, the shallow water wave array data initially processed at 
Scripps Institute of Oceanography, and the processing of the drifter and GPS data 
acquisition are described. 

1. Airborne Radar Data 

The onboard high density digitally recorded radar data is processed into 
a two dimensional imaging map at JPL. The Airborne Imaging Radar (AIR) off-line 
processor is run on a VAX 11/730 minicomputer assisted by an array processor FPS- 
5210. The AIR processor is able to provide two main products: a survey product 
and a high resolution product. The survey processor is a quick look processor, 
scanning the raw data tape to determine the quality of the data and whether the 
data is suitable for futher processing. The high resolution processor provides a high 
quality complex image. The functions which were used to process the data in the 
present study are as follows (Nguyen, 1987, Held et al., 1988). 
a. Range Compression 

The range complex matched filter to the transmitted signal is gen¬ 
erated from a real range offset signal whose spectrum is the complex conjugate of 
the transmitted signal. In order to suppress sidelobes, the range reference function 
is weighted by a cosine function. The range compression is done by inverse Fourier 
transforming the product of two Fourier transform sequences of the match filter and 
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the data. In order to reduce spectrum overlap, the data were subjected to a four-fold 
decimation by applying a bandpass prefilter in the azimuth dimension centered at 
the Doppler centroid frequency. This decimates the azimuth resolution by a factor 
of 4 relative to the nominal SAR resolution. 

b. Azimuth Transform 

The temporal range history, that is, the phase history of a station¬ 
ary synthesized target, is a function of the instantaneous range difference between 
the airborne radar and the target. Therefore, the range of this target can be approx¬ 
imated as a quadratic function of time over the synthetic radar aperture. This effect 
is known as range curvature, or range migration, and has to be corrected so that the 
energy of a point target will be contained only in range bias parallel to the azimuth 
dimension. This correction is done in the frequency domain. The range compressed 
data are first read into azimuth direction format and then Fourier transformed using 
a 8192 point FFT. 

c. Range Cell Migration (RCM) Correction 

The RCM correction algorithm rearranges the data so that the 
quadratic phase history effect in the range dimension is mitigated. 

d. Azimuth Compression 

The data is compressed in the azimuth dimension by convolving 
the data with an azimuth matched filter. This operation is done in the frequency 
domain by taking the FFT of the azimuth matched filter, multiplying it with the 
FFT of the range compressed data, and taking the inverse FFT to produce the one 
look complex image. 

For each of the three legs processed for this study, four separate complex 
image files of the same scene for each antenna are received from the high resolution 
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AIR processor. The size of these files are 750 range lines with 1024 azimuth pixels 
in each line. The initially processed pixel spacing are approximately 6.7 m in the 
viewing direction (range) and 3 m in the azimuth direction. In order to reduce 
speckle noise, the corresponding pixels of each set of the four files are averaged 
incoherently in azimuth. The resulting image from each antenna has therefore a 
degraded azimuth pixel size by a factor of 4, namely 12 m. These two sets of 
files are further registered and combined pixel by corresponding pixel to obtain 
a phase and intensity INSAR image. The phase of each pixel in the combined 
interferogram is the phase difference between the corresponding pixels obtained 
by the two antennas, while the magnitude image represents the product of the 
amplitudes of the corresponding pixels. A brief description for obtaining the INSAR 
images is depicted in Appendix A. 

Further, the aircraft yaw and pitch, which cause sideways motion errors 
were subtracted from the interferogram resulting in pixel phases which represent 
spatial ocean surface velocity fields only. 

The image processed from each flight leg (see Figure 4.5) was received 
from JPL on computer compatible 6250 BPI magnetic tapes. The analysis of 
these images was done at the Air-Ocean Interactive Digital Environmental Analysis 
(IDEA) Lab, which is a remote sensing image processing laboratory for geophysical 
data at the Naval Postgraduate School. 

2. Pressure Array Data 

Time series records of 2048 samples of pressure data from all four chan¬ 
nels (sensors) were processed for the required experimental period 1000-1500 PDT. 
These records were examined for data gaps and signal interference. Some of the 
relevant recognized errors which were taken care of in the initial preprocessing by 
Scripps Institute of Oceanography are: Spikes — Any value of the time series that 
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is more than five times the standard deviation from the mean is considered as a 
spike. Spikes are replaced with the previous value in the time series. Absence of 
zero crossing — a wave time series (after detranding and demeaning) that does not 
cross the zero mean level for a specified number of points is considered unacceptable. 
Filtering — the tidal component is removed. 

The initially processed data by Scripps Institute of Oceanography was 
recorded on a 9 track magnetic tape. Further data analysis was done at the Naval 
Postgraduate School (NPS) and described in the next chapter. 

3. Lagrangian Drifter Data 

Ranges from up to four Mini-Ranger shore stations were logged onboard 
the R/V Ricketts on a mini-computer at 5-second intervals while the vessel ma¬ 
neuvered slowly alongside each drifter. The positioning records indicate that vessel 
speeds while passing the drifters were typically less than 2 m/sec. Hence, vessel po¬ 
sitions were measured at a spacing of 10 meters or less alongtrack (during the Group 
One deployments, data logging was at 10-second intervals which implies positioning 
data about every 20 meters alongtrack). 

Times when each drifter was abeam the Mini-Ranger antenna were hand- 
logged for later interpolation of vessel position from the Mini-Ranger data file. The 
drifters always were passed down the starboard side at a distance of about 2 to 3 
meters from the vessel (antenna). 

Drifter positions were assumed to be the same as the vessel position, 
linearly interpolated between successive Mini-Ranger fixes to the time when the 
drifter was abeam. This procedure does introduce some uncertainty in drifter po¬ 
sition, because (1), the drifter actually was 2 to 3 meters off the antenna position, 
and (2), the interpolation assumes the vessel followed a straight line during the 5 to 
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10 seconds between recorded fixes. However, these assumptions do not practically 
degrade the surface velocity estimates as discussed in the next chapter. 

4. Global Positioning System (GPS) Data 

A set of eleven ground points, visible in the SAR images, were selected 
and positioned to determine their ground coordinates in latitude and longitude. 
Typical points were road intersections and fence line corners. Because the pixel 
size on the images being used was on the order of 20 1 meters, it was felt that the 
positioning technique did not need geodetic survey quality precision. A single GPS 
receiver was Ut>ed with the antenna set as close as practical to the imaged point. 
Once the receiver had acquired the satellites and the readings stabilized (about 15 
minutes), a single position and PDOP value (precision indicator) was recorded and 
the equipment moved to the next point. 

The GPS receiver was operated in the point positioning mode, tracking 
at least four satellites. It should be noted that the observations were made in early 
December 1989, using Block I satellites before “selective availability” was activated, 
which otherwise could have degraded the accuracy of this point positioning technique 
significantly. The standard error in ground coordinates was estimated to be 30 
meters, based on prior experience using GPS for geodetic surveys in the area. 

The ground coordinates obtained from the GPS receiver were latitudes 
and longitudes on the World Geodetic System 1984 (WGS 84) datum, essentially the 
same as the North American Datum of 1983 (NAD 83). Knowledge of the datum 
was important because the bathymetric data was provided on an earlier datum, 

NAD 27, and the shift between the two is significant. 

'The rectified image is of size 1024 x 1024 pixel (see Chapter VII. B) monitored cn a screen of 
512 x 512 (lines x columns). Therefore, the azimuth and range pixel size are doubled (e.g., 12.1 
x 2 m in azimuth) due to this constrain. 
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Elevations of the ground points also were needed, to compute a small 
correction to the range ( y ) pixel numbers. The elevations were scaled from a to¬ 
pographic map of the area, specifically the MARINA quadrangle published by the 
United States Geological Survey. Image coordinates for each ground control point 
were scaled from compressed versions of the three images. It was necessary to use 
compressed images in order to fit the original 1024 x 1024 pixel format on the 
512 x 512 display available. The standard error in picking image coordinates was 
estimated to be three pixels (0.6%) on the 512 pixel format. 

Of the eleven ground control points originally selected, five were dis¬ 
carded because of unresolved blunders, most probably mistaken identification of the 
ground point in the image. This left a set of six ground points to use in the adjust¬ 
ments. A listing of these points and their coordinates is given in Table 1, Appendix 
C, along with other constants pertinent to the adjustment process. 

5. Summary 

This chapter has described the methods, asset performance and data 
collection of the Wave Field and Surface Currents Marina 1989 experiment. The 
digital remote sensing data sets of the L-band INSAR and SAR images from the 
orthogonal and/or opposite legs are further processed to analyze surface currents 
and wavefields. These results will also be compared with conventional data sets 
processed from in situ observations of drifters and wave array data (see Chapter 
VII). 
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VII. DATA ANALYSIS METHODS AND 

RESULTS 


A. INTRODUCTION 

Methods by which the raw data were processed to obtain radar remote sensing 
spatial images and in situ temporal records were described in the previous chapter. 
In this chapter, the SAR and INSAR digital images are further processed to obtain 
two-dimensional wave number spectra of dominant waves, nearshore surface current 
vectors, and estimations of the decorrelation time of the backscattering ocean sur¬ 
face. INSAR wavenumber spectra (wavelength and direction) are compared for the 
first time with in situ wave array spectra and theoretical refraction and shoaling 
predictions. The INSAR measured ocean surface current field is compared with 
Lagrangian drifting buoys following surface currents. The results are in agreement 
and demonstrate the feasibility of INSAR to measure ocean surface wave spectra 
and currents. 

B. SAR AND INSAR IMAGE PROCESSING 

The SAR initial processing at JPL basically generates the two-dimensional 
digital image with a number of inherent geometric and radiometric distortions. Usu¬ 
ally a number of corrections needs to be applied in order to have a correct represen¬ 
tation of the imaged scene. The major post processing methods which were applied 
to the data presented in this study are described below. Digital image data sets 
obtained from flight paths (Figure 6.5) in the Marina experiment September 8,1989 
are presented. 
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1. Geometric Rectification 

Since the image has equal pixel spacing in the slant range, the data 
have to be interpolated in order to get equal pixel spacing in the ground range 
format. The raw image of 750 lines (in range direction) of equal slant spacing and 
1024 pixels in each lines (in azimuth direction) is first scanned to find the number 
of initial unprocessed lines in the image provided by JPL. Next, the meaningful 
nonzero lines are linearly interpolated with equal slant range pixel spacing into 1024 
equal ground range spacing. 


A'to) = A (S<) + 


d(Si+i) - A(Si) 
Si+l — Si 


C 9i - St) 


i = 1,2,. ..,750 
j = 1,2,...,1024 


(7.1) 


where A'(gj) is the new interpolated pixel value, S',-, g, are the slant and ground 
range coordinates of a given pixel with a fixed azimuth coordinate before and after 
rectification, and A(S{) is the original pixel value centered at 5,-. Since the images 
cover relatively small regions (~ 12 km x 7 km) the effect of Doppler shift due to 
the Earth’s rotation which might skew the image in the azimuth dimension is not 
considered. 


2. Radiometric Rectification 

As described in the previous chapter, each image processed data set 
obtained from the three different flight lines in the Marina experiment consists of: 
two complex SAR images and two interferometric 1024 x 750 phase and amplitude 
images. The amplitude or brightness of the INSAR image is similar in a way to the 
conventional SAR image except that it was obtained by taking the modulus of the 
product of corresponding complex pixel images from the two antennas. The image 
pixel brightness /,• is a function of the reflectivity or the backscatter cross section c 
of the corresponding imaged surface area 


u = /(*•), 


(7.2) 
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where /(<r,) is a function of the radar transmitted power, the antenna gain pattern, 
the effect of range attenuation, and the transfer function of the radar receiver and 
the post processor. For an airborne imaging radar (low altitude, small swath) it is 
assumed that the main contribution to brightness attenuation is the range parame¬ 
ter, which attenuates the two way transmitted signal by the fourth power R~ 4 . In 
order to mitigate this effect in the amplitude (or intensity) images, each range line 
is normalized in the image pixel brightness by 

'MD* < 7 - 3 > 

where R\ is the slant range to the first line and 7?,- is the range of line i (t = 
1,2,., 1024). This operation yields a more uniform brightness image in the range 
direction. Figure 7.1 illustrates two brightness images from the northern flight path 
360° before and after radiometric rectification; the difference between the two images 
is clearly observed. Figure 7.2 shows a rectified brightness image from the western 
and eastern flight lines, legs 270° and 90°. 

3. Computation of Image Coordinates 

Ground coordinates for each flight path for SAR, and hence INSAR, 
images were found by computing parameters of a best-fit transformation between 
image and ground coordinates using a small set of GPS surveyed ground points 
visible in each image. A least-squares, minimum-variance adjustment technique 
was used to solve for the five parameters describing a simple rescaling rotation and 
translation transformation: 

k' x = pixel size scale in the azimuth direction (along track); 
k' = pixel size scale in the range direction (cross track); 
a = geographic azimuth of the flight line; 
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Eo = easting ground coordinate of image origin; and 
No = northing ground coordinate of image origin. 

Hence, the image coordinates, expressed as (x,y) pixel number pairs, were treated 

as an orthogonal coordinate system with constant, although different metric scales 

in the x and y directions, the azimuth and range directions, respectively. 

♦ 

A priori values for the first three parameters k' x , fc/ y , a were provided with 
the image data from JPL, or alternatively were computed from information provided. 
These a priori values were used as additional “observations” in the adjustment 
process, and were accepted as the best values for computations using image data. 

Consequently, the adjustment process served two purposes. First, it 
confirmed that the a priori pixel scale sizes and flight path directions (azimuth 
values) were reasonably correct. Second, the adjustment estimated values for the 
ground coordinates of the image or origin, which was needed in order to register each 
image against the bathymetric data set, surface current ground truth measurements, 
ray traces, and station positions (Figure 7.9). 

4. Interferometric Phase Images 

The interferometric phase images (interferograms) are further processed 
to reduce the speckle noise and to obtain image pixels with equal size in azimuth 
and range dimension. The phase image contains the radial component of ocean 
surface velocity. This physical information has to be carefully preserved throughout 
the image; for instance, smoothing methods of the speckle noise have to consider 
this requirement. The extreme dark and bright pixels representing high positive 
and negative velocities are screened out carefully after determining the dynamic 
range of image pixel values. The upper and lower limits of this dynamic range is 
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Figure 7.1: Conventional synthetic aperture radar image of the Marina 
Beach region in Monterey Bay; Top: before radiometric rectification; Bot¬ 
tom: after radiometric rectification. The aircraft attitude was Eastward 
(90°). 


119 


Figure 7.2: Conventional SAR image of similar area as Figure 7.1 with 
aircraft attitude in: Top: Northward (360°) and Bottom : Westward (270°). 
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determined by scanning the dark and bright pixel values over the image relative 
to the “normal” pixel values. Any pixel image value which is beyond or under the 
determined threshold is replaced by a mean pixel value of its four or eight nearest 
neighbors in azimuth and range dimension. The interpolation of pixels for equal 
size in both dimensions assists the extraction of wavenumber in the two-dimensional 
wavenumber spectra (see below). Figures 7.3, 7.4 are phase images corresponding 
to the subswaths illustrated in Figure 6.5. 

C. MEAN SURFACE WIND DIRECTION 
ESTIMATION 

Fortunately, the imaged scene in the Marina experiment included the Salinas 
River. The phase images (Figures 7.3, 7.4) illustrate the difference between the 
stationary coast and moving water surface. This contrast is significantly enhanced 
as compared with the conventional SAR image. The image of the coast is nearly 
uniform with the expected zero velocity except for the moving surface of the Salinas 
River. As can be seen in Figures 7.3, 7.4, and in more detail in Figure 7.5, the river is 
blocked off from the ocean at this time of the year and constitutes a basin of standing 
water in which waves are only generated by local wind with no contamination by 
long fetch dominant waves or advecting currents. The apparent nonzero surface 
velocity in this basin wes measured and compared with Bragg phase velocity and 
wind direction. Anemometer wind measurements at the time of the radar airborne 
overflights taken close to the shoreline of Marina Beach at a height of 27 m show 
clearly an eastward wind component of almost 2 m/sec. This indicates that for leg 
360° the radar was looking at upwind Bragg waves. The measured mean surface 
velocity obtained from the Salinas River is estimated by INSAR phase difference 
expressed by 
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Figure 7.3: Radar interferograms off Marina Beach in Monterey Bay, 
September 8, 1989. Brightness represents surface motions. The aircraft 
flight path is from right to left along the top of the image (Northward 
(360°} attitude). A narrow band long crested swell from south-southwest 
togetner with a wider band wave field from north-northwest are clearly 
observed as they propagate ashore. Topi before filtering; Bottom : after 
Gaussian filtering. 
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Figure 7.4: Radar interferograms of the similar area as in Figure 7.3 
but for flight paths: Top: 90° and Bottom : 270°. 
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Figure 7.5: Salinas River interferograms: Top : for flight path 360°; Bot¬ 
tom: for flight path 270°. 
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2tt B 
X V 


U r 


(7.4) 


The radial velocity U r is in this case directly proportional to the net contribution 
by the Bragg wave phase velocity propagating towards and away from the radar 
in the viewing direction. The radial component of the phase velocity contribution 
is rectified by the local incidence angle to obtain surface components in the range 
direction: 


C P = - 


U t 


sin 7 


(7.5) 


The phase difference a, which can vary between 0 and 2n radians, is repre¬ 
sented by 0-255 image pixel gray shades with zero velocity at pixel value of 128. 
The estimated surface phase velocity CpsR in the Salinas River from the image pixel 
values is given by 

-) TFriizz. (™) 


* /PIXVAL - 128' 
C'Psn — 1 — 


128 J 2 B sin 7 

where PIXVAL is the actual pixel image value. From linear theory, the phase 
velocity of a surface gravity wave is taken to be: 


<Wi 


(7.7) 


Substituting the Bragg resonant condition (Equation 2.6) into Equation (7.6) results 


in 


°p = 


9 X 1 


9 9 • ( 7 - 8 ) 
y Z 7 T 2 Sin 7 

Using Equations 7.5, and 7.7, the estimated apparent INSAR Bragg waves phase ve¬ 
locity is computed from two observations and compared with theoretical prediction. 
The results are shown in Table 7.1. 
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TABLE 7.1: Surface Bragg Wave Phase Velocity Estimated from INSAR 
Interferogram and Linear Theory Phase Velocity 


Flight 




Path 

7 (deg] 

(?psr [cm/sec] 

Cpth [cm/sec] 

360° 

46 

53 ±2 

51 

270° 

33 

—41 1 ± 7 

58 


1 The minus sign indicates Bragg waves propagating away from the aircraft in range direction. 

Based on these results, the surface wind is infered to have had mainly an 
eastward component with a weaker ambiguous northerly component. The result 
clearly supports the assumption that in flight path 360° the radar was looking at 
an upwind surface roughness while in flight path 270° the radar looked at a weakly 
crosswind scene. This indicates that the resonant Bragg waves responsible for the 
radar reflectivity travelled ambiguously in both directions (for flight path 270°), 
towards and away from the radar, resulting in difficulties in estimating their apparent 
contribution to the observed radial velocity from the ocean surface. 

D. CORRELATION OF OCEAN SURFACE 

REFLECTIVITY WITH WIND DIRECTION 

The estimate of the surface wind pattern from the Salinas River interfero- 
grams (Table 7.1), are further infered regarding the mean surface reflectivity ob¬ 
tained from the conventional SAR amplitude images from the three flight paths. 
The surface wind direction as estimated in the previous section had mainly an east¬ 
ward component with a weaker component southward. The image intensity, which 
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is a function of the backscatter cross-section of the corresponding surface, should in 
turn be most intense when looking upwind and weaker when looking crosswind. 

Corresponding mean image intensities from the three flight paths were taken. 
The results are not surprising and are consistent with the estimates of the Bragg 
phase velocity in Salinas River. 


7360° 
•* 90* 

7 90* 

■f270® 


= 4-i-7 


= 1.5-7- 3 


(7.9) 


(7.10) 


These results clearly show that the surface reflectivity is strongest when looking 
upwind, weaker when the radar is looking cross-upwind, and weakest when looking 
cross-downwind. This might be one of the reasons for the difference in the quality 
of the phase images of opposite flights (90° and 270°) covering similar scenes (see 
Figure 7.3). The eastward 90° interferogram seems to be less noisy because the 
essential backscatter signal is sufficient enough to extract the phase information. 
The radar signal obtained from the Bragg backscattering waves is probably less than 
the minimum required in the 270° flight path to obtain a meaningful interferogram 
representing ocean surface dominant waves and currents. 


E. ESTIMATION OF SCENE COHERENCE TIME 

The estimation of the radar scene coherence time of the ocean surface is 
processed following the method described in Chapter IV. D2. Image intensities of 
similar subscenes obtained separately from the two antennas and from the product 
of the two were averaged. The scene coherence time was computed using Equation 
(4.61) 

T c = Jtt (—21n/<y/<o) -1 ^ 2 
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The scene coherence time was computed for all the three flight paths for different 
subscenes. The corresponding averaged subscenes were chosen from near to far 
incidence angle. The results are depicted in Table 7.2. For the Northern flight path 
(360°), the viewing ocean scene is only between 40°-50° of incidence angle; the scene 
coherence time first increases slightly (until 46°) and then remains almost constant. 
For further discussion of these results, see Chapter VIII. B. 

F. AMBIENT SURFACE CURRENT ESTIMATES 
FOLLOWING DRIFTING BUOYS 

Vessel positions were computed from the Mini-Ranger ranges using a least- 
squares technique called “variation of coordinates”. Descriptions of the procedure 
can be found in most surveying texts. The presentation in Cross (1981) is virtually 
identical to that used in the present analysis. 

Although the Mini-Ranger system computes and logs a position on line, all 
range data were reprocessed in order to eliminate errors caused by erratic reception 
from the station at MBARI RM (Table 6.6). In fact all data from that station were 
eventually rejected. 

1. Positional Accuracy 

The two-dimensional random errors in vessel position are estimated as 
part of the “variation of coordinates” procedure. These errors are primarily func¬ 
tions of the relative geometry between the ship and shore stations, and the esti¬ 
mated precision of a single range measurement. This random error is expressed as 
a variance-covariance matrix computed as follows: 

C = a 2 0 {A T PA)~ 1 (7.11) 
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Table 7.2: Estimates of Scene Coherence Time as a Function 
of Incidence Angle Obtained from Three Flight Paths in the Marina 
Experiment, September 8, 1989 


Mean Incidence 

Scene Coherence Time r [msec] 

Angle 7 

leg 90° 

leg 270° 

leg 360° 

25° 

136 

121 

— 

29° 

146 

123 

— 

33° 

139 

134 

— 

36° 

117 

116 

— 

o 

O 

121 

105 

119 

42° 

114 

97 

134 

45° 

i 

96 

95 

133 

o 

CO 

^5* 

109 

95 

138 

50° 

112 

94 

140 

Mean 

121 

109 

132.8 

Standard Deviation, a 

15.2 

13.2 

7.35 


129 




where C is the 2x2 symmetric variance-covariance matrix with the diagonal elements 
being the variances of easting and northing coordinates and the ofT-diagonal being 
their covariance. The quantity Vq is the “a priori variance of unit weight” which 
scales C. Practically, trjj is assumed to be unity in these computations and checked 
with an “o posteriori 1 '' estimate derived from the actual range measurements in each 
fix. 

The matrix A, is the N x 2 Jacobian for the transformation from the 
iV-dimensional range coordinate system to the two-dimensional UTM coordinate 
system. In this case N typically was 3 for the ranges from reference stations WATS, 
HAYS, and DOP RMl finally used (see Table 4.6). 

Aij = {Xj - Sij)/Ri t = 1,2,3; j = 1,2; (7.12) 

where X is the ship position, Si is the z-th station position, and Ri is the horizontal 
distance on the UTM grid between X and Si, 

The matrix P is the N x N weight matrix which is the inverse of the 
variance-covariance matrix for the N range observations. It is assumed the range 
measurements are all independent with equal variances (it is assumed here o = 3 
meter., for any single range measurement). Hence P becomes a simple identity 
matrix times the scalar quantity 1 /a 2 , or 1/9. 

Because the drifters were ofTset from the Mini-Ranger antenna, an ad¬ 
ditional random error must be included. The ofTset distance was between 2 and 3 
meters, but the direction was not known. In theory, the course-made-good derived 
from successive fixes could be used to estimate the ship’s heading, and hence the 
direction to the drifter. However, uncertainty in direction still would be very large 
because of the relatively slow speeds and inherent “crabing”. It was simpler to just 
add a circular random error of 2 meters to the variance-covariance matricies derived 
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above. It can be shown this is accomplished by simply increasing the easting and 
northing variances by 4 m 2 each. 

This additional circular random error also serves to represent an uncer¬ 
tainty due to the linear interpolation assumption. No theoretical approaches could 
be found in estimating the random error in the interpolation, however, based on 
simulation runs, the maximum error due to interpolation between fixes of 10 m 
apart would be about 2 m. 

In Appendix D, the variance-covariance matrices computed for each ob¬ 
servation of drifters from Group One and Two are presented. Only one matrix was 
computed for each group of observation because the drifters were so closely bunched 
that the relative geometry did not change significantly and hence, the random errors 
were essentially the same for all. 

In general, errors increased slightly closer towards the shore. Circular 
errors increased from about 11 m to 13 m (CE 95%), going from the most offshore 
positions where Group One was deployed to the nearshore area where Group Two 
was recovered. Also note that the northing coordinate was determined more pre¬ 
cisely than the easting, by a factor of about 3. Both effects are due to the fact that 
the reference shore stations lay essentially along a north-south line passing close to 
the experiment region. 

2. Accuracy of Velocity Estimates 

Drifter velocities were computed by the difference of positions divided 
by time interval. Standard error propagation techniques were used to estimate a 
variance-covariance matrix for the U and V velocity components. The reader is 
refered to the MATHCAD documentation for a description of the math and the 
actual computation of the estimates given below. 
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The error estimates of the sorted data are small, of the order of 0.1 
cm/sec as shown in Table 7.3. The errors depicted are small primarily because 
of the relatively long time intervals between drifter position observations, typically 
about 40 minutes. The relatively large time period of drifter positioning is allowable 
since the environmental conditions on September 8,1989 were almost steady for this 
time interval. Also note that these error estimates assume that U and V components 
are constant over the period between drifter positionings. In other words, these are 
estimates of the error in the mean velocities. 


TABLE 7.3: Standard Deviation of Velocity Components (cm/sec) 




°U 

crv 

Group One: 

Obs 1 to 2 

0.13 

0.08 


Obs 2 to 3 

0.19 

0.11 


Obs 3 to 4 

0.22 

0.12 

Group Two: 

Obs 1 to 2 

0.14 

0.07 


Obs 2 to 3 

0.20 

0.10 


G. INSAR MEASUREMENT AND GROUND 

TRUTHING OF NEARSHORE OCEAN SURFACE 
CURRENTS 

INSAR phase images (Figures 7.3, 7.4) were further processed to obtain esti¬ 
mates of very near ocean surface currents. The radial component of surface velocity 
U x of each pixel is rectified from the incidence angle dependence resulting in the hor¬ 
izontal component of ocean surface velocity in the viewing (range) direction. The 
bias introduced by the Bragg wave phase velocity estimated from the Salinas River 
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(see Chapter VII. C) is then subtracted. The resulting phase image is averaged 
two-dimensionally to mitigate the quasi-periodic wave orbital velocity. The scale of 
averaging is the orbital motion scale of the largest observed wave. All smaller wave 
multiples of this scale are averaged out. The surface velocity image from one flight 
path represents only one component of the surface current vector in the viewing di¬ 
rection. Two orthogonal flight paths observing the same area are needed to estimate 
the full current vector. The phase image data from flight paths 360° and 90° were 
selected to obtain interferometric virtual estimates of ocean surface current fields. 
These estimates are compared with simultaneous observations of the Langrangian 
floating drifters, and good agreement is obtained. 

1. Additional Considerations and Processing 

When computing current estimates from INSAR phase images, the fol¬ 
lowing additional considerations were taken: 

• When using the interferometric SAR technique to measure surface velocities, 
one has to consider sideways motion of the antennas, introducing a phase shift 
making the imaged scene appear to move when it is stationary. The prominent 
errors in the image phase data were corrected at JPL (see Chapter VI. F) using 
the Inertial Navigation System (INS) data and land imagery as a stationary 
reference (Goldstein et al., 1989). 

• Mini-Ranger positions of the floating drifters were transformed to image pixel 
coordinates using the GPS adjustment method (see section B3). Image sub- 
areas from flight paths 90° and 360° containing the observed path of the La- 
grangian drifters were selected to match the current estimates as spatially 
similar in nature as possible to the in situ observation. 
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• The image constructed from the 360° flight path looking upwind is less noisy 
compared with images constructed from flight paths 90° and 270° looking 
crosswinds (see section C, D). Therefore, prior to the computation of current 
estimates, the phase image from the eastward flight path was smoothed. Pixel 
phase values, which deviate by more than two standard deviations from the 
mean phase value of a selected subarea, were rejected, resulting in a less noisy 
phase image for current computation. 

An additional problem is the appearance of slicks (seen as dark regions 
in the conventional SAR amplitude image, Figures 7.1, 7.2) in the selected area. In 
these areas the essential backscattering amplitude is probably too weak, i.e., very 
noisy, to extract a reliable interferometric phase due to surface motion. Nearby 
regions with no slicks were selected to compute current estimates for areas observed 
by drifters 3, 4. 

2. Comparison of INSAR Current Estimates with Lagrangian 

Drifter Measurements 

The simultaneous in situ surface current estimates are shown in Figure 
7.6. Three groups of five Lagrangian drifters each were deployed close to the time of 
INSAR aircraft overflights (see Chapter VI. C). Radar data was initially acquired 30 
minutes before the predicted time of maximum flood tide current (1320 PDT). This 
timing was predicted on observing the stronger flows to make the current estimates 
more reliable in an area (Marina Beach) which is known to have unpredictable weak 
tidal currents. 

The Lagrangian drifters (Figure 7.6) experienced a dominant eastward 
component and a much weaker northward component. Some estimates of surface 
current derived from interferometric SAR data and Lagrangian drifters are tabulated 
in Table 7.4. It is shown that the dominant eastward components are in better agree- 
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LONG 

Figure 7.6: Direct observation of surface currents following Lagrangian 
drifters near Marina Beach in Monterey Bay, California. Th*. arrows 
point in the direction of the current and the length of each arrow is 
proportional to the observed t 'eed. 

ment (correlation 0.9) than the weaker and hence, difficult to determine, northward 
component (correlation 0.5). The velocity magnitudes are most consistent with a 
correlation of 0.94. The reasons for these results are assumed to be associated with 
the physical difference between the two types of measurements, Lagrangian versus 
spatial synoptic observation, and are further discussed in Chapter VIII. C. 

The drifter velocity estimates with almost nil standard deviations (a u = 
0.18 cm/s, c T — 0.1 cm/s) were discussed in the last section. The INSAR mea¬ 
surements had a larger standard deviation of 20-30 cm/sec computed from selected 
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TABLE 7.4: INSAR Estimates and In Situ Observation of Ocean Surface 
Currents (in cm/sec). Subscript d Denotes Drifter and s for INSAR. 


Drifter 

Zonal 

Component 

Meridional 

Component 

Velocity 

Magnitude 

Number 

Ud 

u, (360°) 

Vd 

w. (90°) 

v d 

u a 

1,2 

20 


5 


21 

m 

3,4 

18 

15 

5 

4 

19 

11 

5 

18 

16 

4 

5 

18 

Efl 

6,7 

22 

17 

5 

9 

23 

19 

8 

21 

17 

5 

8 

22 

19 

9,10 

22 

18 

4 

5 

22 

19 

Average 

20 

17 

5 

6 

21 

18 

rms error (cm/sec) 


0.5. 


1.7 


0.5 


subareas, which were averaged to estimate INSAR current components. This large 
velocity deviation is believed to be primarily associated with the orbital velocity of 
the quasi-periodic surface gravity waves which are averaged out. To get an estimate 
of the order of magnitude of this assumption, consider a a deep water surface gravity 
wave propagating in the cross-track (range) direction of a moving radar platform. 
For this geometry, the observed radial component is the actual wave orbital veloc¬ 
ity. For waves propagating in the along track (azimuth) direction, only the vertical 
particle velocity is observed. Therefore, the largest observed orbital velocity in deep 
water is given by 


U _ » 9T' 

VoB ~lT 


(7.13) 


where II is the wave height. The orbital velocity of a shoaling wave increases with 
decreasing wavelength; hence, the expected dominant contribution to the observed 
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mean orbital velocity is from the shorter observed swell component (see Chapter IV. 
C). 

For a wave with L = 97 m, period V — 9.1 sec and a rms wave height 
H a = 0.32 m (the wave height and period are estimated from the in situ power 
spectrum, see Figure 7.10) as observed in the Marina experiment, U t ob = 0.20 
m/sec. This result is consistent with the computed standard deviation from the 
estimated mean surface currents. It indicates that most of the contribution to the 
standard deviation is from the quasi-periodic orbital velocity of dominant waves 
which are averaged out. Additional contribution to the standard deviation of the 
mean surface current component estimates is attributed to speckle noise. Estimates 
of current velocity vectors for radar and direct observed measurements are shown 
in Figure 7.7. The agreement is felt to be encouraging. For further discussion on 
these results, see Chapter VIII. C. 



50 CM/S 

Figure 7.7: Detailed presentation of Figure 7.6 together with INSAR 
current estimates (dashed arrows). 






H. OBSERVATION OF SPATIAL AND TEMPORAL 
NEARSHORE OCEAN WAVE SPECTRA 

Shoaling and refraction of waves can be characterized by two dimensional 
spectral analysis. The two-dimensional wavenumber (k = ( k xi k v )) spectra obtained 
by INSAR imagery are tested against the conventional measurement of directional 
spectra computed from time series records obtained from the pressure wave array 
in the shallow water off Marina Beach. In addition, the two swell components 
observed by the shallow water array were used as inputs to compute wave ray traces 
using Dobson’s numerical model (see Chapter V. B). Two-dimensional INSAR wave 
number spectra were computed for subscenes along the ray traces. The results 
are consistent with the refraction model outputs and with the shallow water array 
directional spectra. 

1. INSAR Two-Dimensional Wavenumber Spectra 

Following the description in Chapter V. C, a software package was de¬ 
veloped to calculate two-dimensional wavenumber spectra from INSAR images. 

The rectified phase images illustrated in Figures 7.3, 7.4 were further 
processed by linear interpolation (see Equation 7.1) to equalize the azimuthal pixel 
size with the range size resulting in square pixels. This operation was done to 
simplify the computation of the wave lengths and directions of the peaks of the 
wave number spectra. 

Due to the complex shallow water region in which the experiment was 
conducted, relatively small subscenes had to be selected to fulfill the homogeneity 
requirement for spectral analysis. Trade off is between satisfying the requirement 
of homogeneity and obtaining confidence in the spectral estimates. To obtain con¬ 
fidence in the spectral estimates, a minimum sea surface area equivalent in size to 
approximately 6-10 wavelengths on a side must be used to calculate INSAR image 


138 







spectra to provide reasonable estimates of wavelength and direction. Therefore, a 
flexible 2D FFT code was developed which allows a selection of any number of pix¬ 
els in either dimension resulting in rectangle sub-areas for spectral analysis of any 
size. If the number of pixels in the sub-area is not a power of 2, the sub-area is 
zero paded to the nearest power of 2 number of pixels. This increases artificially 
the effective resolution for wavelength and angle computations. The edges of the 
subscenes were cosine tapered to reduce spectral leakage and demeaned. Figure 7.8 
shows some of the sub-areas that were subjected to spectral analysis. The INSAR 
image spectra of the corresponding sub-areas are illustrated in Figures 7.12 to 7.15. 
They are smoothed with a Gaussian smoothing kernel providing spectral estimates 
with 25-30 degrees of freedom. 

2. Shallow Water Array Directional Spectra 

Computation of shallow water array directional spectra was done us¬ 
ing the Grauzinis (1989) method (see Chapter V. C). The available Fourier coeffi¬ 
cients of the wave field angular distribution function were computed using Equations 
(5.31)—(5.35). The model distribution function that was used to match the Fourier 
coefficients is the wrapped Cauchy distribution 


if 1 - m\ 

2tt 11 — 2mi cos 0 + m\ 

This distribution gives suitable narrow estimates and allows for bimodal directional 
spectra. 

The spectral analysis was performed for a one hour time series 1200— 
1308 PDT on September 8, 1989 with subsequent FFT’s computed over 128 second 
intervals with a sampling rate of 0.5 Hz. The averaged spectral estimates have ap¬ 
proximately 64 degrees of freedom. As mentioned in Chapter VI. F2, the data were 
detrended and demeaned before being Fourier transformed to eliminate contamina- 
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Figure 7.8: Radar interferogram of Marina Beach area, September 8, 
1989. The indicated sub-areas were subjected to spectral analysis with 
results appearing in Figures 7.12-7.15. 


tion of the wave frequency bands due to leakage of low frequency signals. The wave 
pressure power spectrum is calculated and transformed into surface elevation and 
surface slopes spectra by applying linear wave theory transfer functions as described 
in Chapter V. C, and subsequently converted from a frequency angle presentation to 
a directional wavenumber presentation via intermediate water linear wave equations. 

3. Comparison of INSAR 2D Wavenumber Spectra with Temporal 
Wave Array Spectra 

Wave refraction rays for the two dominant waves observed by the Ma¬ 
rina shallow water array were constructed using the Dobson model. The rays at 
these peak frequencies (0.063, 0.109 Hz) were propagated offshore. The predicted 
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121° 4 9.6’W 


Figure 7.9: Observed ray traces — from INSAR spectral analysis versus 
predicted wave orthogonals for the two swells measured by the shallow 
water array of Marina Beach. The numbers are the centers (stations) of 
the 2D FFT sub-areas. 
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Table 7.5: Actual INSAR 2D Wavenumber Wave Array Spectra Results and Predicted Refrac¬ 
tion and Shoaling Data 
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wavelength and direction along the refracted ray (about every 800 m) were compared 
with the actual wavelength and direction of the two dimensional INSAR image 
wavenumber spectra from the Northern (360°) flight path image (see Figure 7.8). 
The analysis is summarized in Figure 7.9 and Table 7.5. 

The effects of tides and coastal currents on wavelength and propagation 
direction were ignored and considered negligible when constructing the ray traces. 
The comparison between the predicted and INSAR results clearly demonstrate the 
feasibility of INSAR in imaging ocean wave spectra in a complex region under the 
given experimental environmental conditions. 

An inherent limitation in applying the refraction model to INSAR’s data 
is that the bathymetry used in the model is not as accurate as would be desired. 
The density of depth data varied based on various ship tracks and was available 
only every 6 seconds. 

The in situ wave array spectral results are more reliable for comparison 
with INSAR’s measurements. The surface elevation power spectrum obtained by the 
wave sensors is illustrated in Figure 7.10. The two observed peaks in the spectrum 
correspond to wave periods of T\ = 15.9 sec and T 2 = 9.1 sec. The wavelengths 
(Station 1, Table 7.5) are calculated at the appropriate sensor depth (16 m) using 
the linear wave dispersion relation (Equation 5.3). The higher energy region around 
the short wave peak (T\) shows a range of spectral average (0.1 to 0.17 Hz), while 
the region around the lower energy peak (T 2 long swell) is much narrower; thus, the 
spectrum indicates the existence of two swells with the shorter wave field having a 
broader spectrum. 

The directional wave spectrum computed using the exact Fourier co¬ 
efficient representation method (Grauzinis method) is presented in Figure 7.11. 
The shorter wave field propagates from west-northwest to east-southeast, while 
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the longer swell propagates from south-westsouth to east-northeast. The wave¬ 
length estimates from the spectral analysis of the two data are identical while the 
directions are in good agreement (see Table 7.5). 



Figure 7.10: In situ surface elevation power spectrum measured by the 
Marina shallow water array on September 8, 1989, 1200-1308 PDT, at a 
depth of 16 m. 

In order to reduce bathymetric constraints, two typical sizes of image 
sub-areas were selected for computing the 2D wavenumber spectra. For the shorter 
wave field, sub-areas of ~ 1400 x 700 m were selected for computing wavelength 
and direction depicted in Table 7.5, while for the long crested swell sub-areas of 
~ 1400 x 1400 m with overlap between the sub-areas were used. The variation of 
bottom topography over the subscene changes the wavelength and direction. This 
distortion is reflected in the 2D wavenumber spectra. Table 7.5 illustrates the 2D 
wavenumber spectra of INSAR along the ray traces depicted in Figure 7.9. 

The compu'ed spectra (Figures 7.12-7.15) illustrates shoaling and re¬ 
fraction of the two dominant swells over a wide imaged region (~ 12 x 0 km). For 
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region A, B where the ocean is relatively deeper, shoaling and refraction are less 
pronounced compared to the southern regions C, D. 

Ill spite of the limitations mentioned above, INSAR image spectra veri¬ 
fied by wave array measurement and compared with model predictions are encour¬ 
aging, and demonstrate the feasibility of this unique technique to estimate ocean 
wave spectra even in a complex occonographic region. For further discussion, see 
Chapter VIII. D. 
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Sept 8, 1909 1200-1308 POT 




Figure 7.11: In situ directional wave spectrum taken at Marina Beach 
for the same time interval as Figure 7.8; Top : 3D representation; Bottom: 
contour presentation. 
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Figure 7.12: Top : Shallow and Bottom: deeper water INSAR image 
wavenumber spectra of sub-area Al, A2 indicated in Figure 7.8. The 
circles represent wavelengths as depicted in Figure 7.11. 
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Figure 7.14; Same as Figure 7.13, but for sub-area Cl, C2. The top 

image spectra is obtained from a sub-area close to the Marina shallow 
water array. 












VIII. RESULTS AND DISCUSSION 


A. INTRODUCTION 

The aim of the experimental part of this study has been to verify the capabil¬ 
ities of INSAR to image near surface ocean dynamics. The method was to compare 
ocean wave spectra and surface current measurement obtained remotely by INSAR 
with conventional in situ observations. 

A discussion of the processed results from Chapter VII is presented in this 
chapter. It is shown that the interferometric technique appears to offer considerably 
improved imagery of ocean surface waves as compared to conventional SAR. The 
additional capability relative to SAR, in imaging ocean surface currents, seems to 
be reliable for oceanographic studies. 

B. TEMPORAL RADAR COHERENCY OF THE 
OCEAN SURFACE 

The fluctuations in the short wave reflectivity of the ocean surface during a 
given observation time caases the temporal correlation function of its reflectivity to 
have a width r c . The azimuth resolution p az is degraded as compared to that of a 
stationary scene (see Equation 3.2) to 



In addition to the degradation in azimuthal resolution due to finite scene 
coherence time, the SAR (and hence INSAR) azimuthal resolution is smeared also 
by dominant wave motions. The azimuthal wave pattern translation and the range 
component of the wave acceleration during the integration time primarily degrades 
the azimuthal resolution. For an imaged swell, the wave pattern translation smearing 
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can be corrected by focus adjustment of the processing matched filter (Raney, 1980) 
and therefore will not be discussed further. As was noted in Chapter IV. C, the 
finite scene coherence time determines the effective integration time and hence the 
azimuthal resolution. The L-band radar scene coherence time estimation from image 
intensities (Chapter VII. E, Table 7.2) is of the order of O(10 2 ) msec. 

The range component of acceleration of a deep water gravity wave is given 

by 

a(x,<) = 2 H (jfij ex P p^~j sin(k • x — wt) cos 7 (8.2) 

where H is the wave height, T' is the wave period, and z is the water depth. At the 
surface 2 = 0. The degradation in azimuth resolution due to wave acceleration is 
taken to be (Alpers and Bruening, 1986) 

Paz = | 7 °r(x, i)T=^y 2 H (jp) sin(k -x-ut)T cos 7 (8.3) 

For typical JPL/DC -8 parameters, wave height of 0.7 m, wavelength of 100 m, 
effective integration time of 0.1 sec. and incidence angle of 30°, p az = 0.7 m. 
Therefore, the degradation in azimuth resolution due to long wave acceleration is 
neglected. 

It is important to note once again that the loss in azimuthal resolution due 
to finite scene coherence time associated with the random motion of scattering el¬ 
ements within the resolution cell, results in a consistent azimuth image smear and 
not defocus, and is therefore not recoverable by focus adjustment or integration 
time settings (Lyzenga and Schuchman, 1983). Therefore, the estimation of scene 
coherence time remains a key parameter in determining the azimuthal cutoff wave¬ 
length of SAR and INSAR azimuthal resolution. If the coherence time r c is much 
smaller than the post processing integration time T, the azimuthal resolution can 
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be approximated (Lyzenga and Shuchman, 1983) to be: 


Paz ~ 


A R 
2 Vr c 


(8.4) 


Based on the results depicted in Table 7.2 and Equation (8.4), the degradation 
in azimuth resolution due to finite scene coherence time for the three flight path 
images in the Marina experiment is illustrated in Figure 8.1. 



Figure 8.1: Degradation in azimuth resolution for the NASA/JPL DC-8 
L-band SAR for the three flight path images in the Marina experiment, 
September 8, 1989. 


These estimates of L-band scene coherence times are consistent with con¬ 
ventional radar observations (100 msec) during the TOWARD experiment (Plant 
and Keller, 1984). The theoretical estimates of the scene coherence time using the 
Pierson-Moskowitz spectrum (Table 4.1) gives larger estimates (3-4 times) com¬ 
pared with the more realistic direct observation estimates. The estimates of the 
ratio K/Kq in Equation (4.61) for all three flight path images were close to 0.9 and 
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larger. This result suggests that for the estimation of the scene coherence time only, 
the measurement sensitivity can be increased by increasing the effective baseline 
separation (9.9 m) between the antennas. Variation of the effective baseline sepa¬ 
ration will allow to investigate the distribution of the temporal correlation function 
which was assumed to be Gaussian in the present study. 

For rougher sea state, it is expected that the observed scene coherence time 
will be shorter. This is because the waves become steeper with higher orbital veloc¬ 
ities resulting in a larger velocity spread of the scattering elements. When break¬ 
ing occurs, the backscattering mechanism is specular in nature (not Bragg). The 
backscattering slopes are propagating at the dominant wave phase velocity (Lyzenga 
et al., 1983). The transformation of the advecting velocity of the backscattering sur¬ 
faces, from orbital to phase velocity when breaking occurs, further reduces the scene 
coherence time. A cusped region is the surf zone, the estimated scene coherence time 
from this region is 30% less relative to estimates from nearby sea surface regions 
with similar incidence angle. 

The difference in azimuth resolution between three flight path images is not 
fully understood at the moment and is probably determined by the speckle and 
random noise and look direction (upwind crosswind). Based on the estimates of 
azimuth resolution (40-70 m) in Figure 8.1, the largest unresolvable azimuth waves 
are of wavelength between 80-140 m (increasing with incidence angle). The orbital 
velocity of these waves modulating the fundamental scattering elements mainly de¬ 
termine the finite scene coherence time (see Chapter IV. C). 

The degradation in azimuth resolution is one of the major reasons for the less 
pronounced wave patterns in the INSAR images obtained from flight paths 270° and 
90°. The azimuthal orientation of the dominant swells makes them less resolvable 
with increasing range as the azimuthal resolution decreases (Figure 8.1 and Chapter 
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VIII. D. 2c). The wave patterns in the 90° and 270° SAR images were more affected 
by velocity bunching compared to the 360° SAR image. This might contribu. j to 
additional azimuthal degradation of the INSAR image. 

In shallow water shoaling and refraction probably also play a role in de¬ 
termining the scene coherence time. Due to refraction, waves tend to propagate 
perpendicular to the bathymetric contours (refraction) which are mainly aligned 
alongshore, while shoaling results in increased wave amplitude and elliptic particle 
orbits. The change in the radar observed orbital velocity due to shoaling is believed 
to be small. The more pronounced effect is probably the refraction effect in which 
wave crests tend to align with the bathymetric contours in shallow water before 
breaking occurs. In deeper water, the wave length and direction are less affected by 
the bathymetry and the direction of wave propagation might be affected by local 
winds. In the Marina experiment the mild wind (< 2 m/sec) had no influence on 
the direction of propagation of the larger unresolved waves (80-140 m) which are 
of our interest since these waves determine the azimuth resolution. Therefore, the 
direction of the larger unresolved waves in deeper water was unknown. For flight 
path 360°, the refraction effect caused waves to be range oriented near the shore 
(smaller incidence angles) and hence the full magnitude of the orbital velocity should 
be observed by the radar, resulting in a faster decorrelation time of the surface. An 
indication of this effect can be noticed through the results depicted in Table 7.2 
for the 360° image. Additional evidence using the present technique and in situ 
observations is needed to get a better understanding of the determination of scene 
coherence time in shallow water. 

An important application of the knowledge of the scene coherence time is in 
defining the optimum physical baseline separation between the two antennas which 
determines the sensitivity of interferometric measurements of ocean surface velocity 
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fields. For a mean estimate of 0.1 sec. of the scene coherence time, it is required 
that the platform travel time through the effective baseline separation (9.9 m) will 
be less than 50 msec. Using the expression for the time delay for interferometric 
measurement A t - r c /2 and L-band NASA/DC-8 parameters yields a maximum 
baseline separation of: 

B m&x = At-2V = 21.4 m (8.5) 

This result shows that the actual baseline separation (19.8 m) used in the Marina 
experiment was almost optimal for the given environmental conditions. For rougher 
sea state and/or higher radar frequencies, the scene coherence time will be shorter 
and hence the baseline separation between the antennas should be reduced. The 
requirement for At = r c /2 probably could be more flexible allowing At < r c and 
hence B mAX could be increased resulting in an improved sensitivity in the INSAR 
phase measurement (given r c > 0.1 sec.). The effects of increasing of the antenna 
baseline separation should be inferred in future experiments. 

C. INTERFEROMETRIC SAR REMOTE SENSING OF 
OCEAN SURFACE CURRENTS 

A primary objective of the present study is the verification of the feasibility 
of INSAR to estimate reliably ocean surface currents with high resolution over large 
regions. The Marina experiment demonstrates this feasibility by comparing INSAR 
estimates with almost simultaneous in situ Lagrangian drifter estimates (see Chapter 
VII. G). A scatter diagram of current magnitude estimates derived from the two 
sensors is illustrated in Figure 8.2. Simple statistical analysis shows a significant 
degree of linear correlation (R = 0.94), with a small consistent bias of 2.5 cm/sec. 
However, due to the relatively small number of compared estimates the statistical 
significance is rather low. 
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DRIFTER CURRENT ESTIMATE (CM/S) 

Figure 8.2: Lagrangian drifter versus INSAR surface currents estimates. 

1. Limitations of Measurement Methods 

In the eastward (90°) flight path image (from which the meridional cur¬ 
rent component was estimated) the radar was looking crosswinds. The estimation of 
the ambiguous phase velocity associated with the Bragg waves was derived from the 
western (270°) flight path image (see Chapter VII. C, G), which included the Sali¬ 
nas River. Unfortunately, the less noisy 90° phase image did not include a nearby 
calm region like the Salinas River to estimate the phase velocity bias. The 270° 
image was not selected for computation of meridional current estimates because it 
was the most noisy image obtained in the Marina experiment (see Chapter VII. D). 
Therefore, the substruction of the velocity bias introduced by the Bragg wave phase 
velocity for the meridional component was not straightforward, resulting in a less 
accurate estimation of this component. 
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Another difficulty that was addressed in estimating current components 
from INSAR phase images was the presence of slicks, which are clearly seen in the 
amplitude S AR images (see Figures 7.1,7.2). These slicks are believed to be observed 
due to local dampening (straining) of the wind generated surface ripples which are 
responsible for the essential Bragg radar backscatter. Ocean surface straining can 
be introduced by various oceanographic phenomenas, one probable driving force to 
create slicks is that the surface signature is related to the convergence of surface 
particle velocities associated with internal waves. The internal waves might be 
generated by tidal flow variabilities over topographic features causing fluctuations 
in the thermocline. The local convergence zones introduce concentration of surface 
organics which dampen the surface ripples due to increased surface tension. Another 
mechanism creating regions of surface convergence are Langmuir cells in which the 
slicks are oriented in the direction of the wind. Because some of the slicks are aligned 
along the bathymetric contours and the wind was weak, it is believed that they are 
related to internal waves rather than Langmuir cells. 

Mean temporal Lagrangian estimates are physically different in character 
from the synoptic instantaneous spatial observation of the surface currents, which 
could be described as a grid of Eulerian snapshots. The Lagrangian drifter trajectory 
would not necessarily follow the spatial current pattern measured instantaneously 
by INSAR, which is evident in Figure 8.2. Averaging of imaged subareas make the 
two different current estimates more comparable. This difference might be one of 
the reasons for the better agreement in the computed current magnitude estimates 
rather than current direction (see Table 7.4). 

2. Current Measurement Conclusions 

The results of the experiment in estimating surface currents with INSAR 
shown in Table 7.4, Figures 7.7, and 8.2, together with the discussion above, suggest 
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a number of conclusions which are basically consistent with those pointed out by 
Goldstein et al., 1989. 

• The method seems to be most reliable when the radar observation is along the 
direction of the strongest current component (flight path 360°). 

• The magnitude, rather than the direction of current estimates, is more con¬ 
sistent with conventional observation, probably because of the fundamental 
difference between the two measurement techniques. 

• The interferometric technique seems to be more accurate when the radar is 
looking upwind (and probably also downwind). In these cases, the Bragg 
waves phase velocity responsible for the radar signal are less ambiguous, and 
the returned amplitude signal is strongest resulting in a better signal over 
noise. 

• The estimates of the bias phase velocity (when looking upwind) are consistent 
with those expected theoretically. 

• The interferogram contains numerous additional features in the ocean flow 
field which were neglected in the present study. Spatial nearshore current 
patterns is one phenomena which should be further analyzed. 

D. OCEAN WAVE SPECTRA 

The results of INSAR image spectral analysis tested against conventional 
wave array spectra estimates and refraction model predictions (Chapter VII. H) 
demonstrate that interferometric SAR under certain geometric orientations appears 
to offer considerably improved (relative to conventional SAR) imaging of waves also 
in complex nearshore regions. 
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1. Results Comparison 

The predicted (Dobson’s model) ray traces expressing the relationship 
between the water depth h and wavelength L find between the depth and direction 
of arrival a (Figures 8.3 , 8.4) are such that the refraction model becomes insensitive 
to changes in h when 

h > 0.4I o , (8.6) 

where L 0 is the deep water wavelength. For the two observed swells in the Marina 
experiment this requirement results in a depth of 54 m for the shorter swell (9.1 
sec) and 163 m for the longer swell. The western edges of the Northern flight path 
(360°) INSAR image which was mainly subjected to spectral analysis (Figure 7.8) 
represents a mean depth of 53 m. This indicates that both waves experience shoaling 
and refraction through the entire imaged scene. 

The constructed wave refraction diagram (Figure 7.9) is consistent with 
the INSAR derived orthogonals. As mentioned previously, wave propagation an¬ 
gles can be sensitive to other oceanographic phenomena such as tides and coastal 
currents. Wavelengths are also subject to these effects but to a much lesser extent 
(Shuchman and Kasischke, 1981, Hayes, 1980). These effects were not considered 
in our refraction model. The observed weak currents (Chapter VII. G) and the 
consistent INSAR spectral results consistent with model predictions, indicate that 
these effects were minor during our experiment. 

Figure 8.5 presents a plot of wave directions Oj, ect determined from 
INSAR 2D FFT versus oi, a 2 calculated from Dobson’s refraction model. The 
plot shows a significant degree of linear correlation (0.88) for the long crested swell 
which experience more significant refraction compared to the shorter observed swell 


160 


o 








I- 1 - 1 - 1 -r - 

220 234 240 240 252 

WAVE DIRECTION L2 (M) 


Figure 8.4: Relationship between water depth h and predicted direction 
of wave arrival or in Monterey Bay based on Dobson’s refraction model. 
Top: for the shorter swell L 0 = 129 m; Bottom : for the longer swell L 0 — 
395 m. 
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(correlations of 0.43). The reason for the better agreement of the narrow band long 
crested swell is probably because this swell is more monochromatic in character 
relative to the shorter swell and hence more comparable to the theoretical model, 
which is based on the analysis of monochromatic waves as input 

The results for INSAR observed wavelengths compared with model pre¬ 
dictions are depicted in Figure 8.6. The linear correlation is better than 0.9 for both 
swells. This is mainly because shoaling is less sensitive than refraction to oceano¬ 
graphic constrains and remote sensing imaging geometry. The measured wavelength 
from INSAR’s image spectra compared with the in situ wave array wavelength esti¬ 
mates and model predictions, is in better agreement than wave direction (see Table 
8.1). 

It is felt that more advanced spectral analysis techniques like maximum 
entropy might improve the data analysis results in shallow water. This and other 
techniques offer the potential of reducing the number of wavelengths required in a 
sub-area which is subjected to spectral analysis (Schuchman and Kasischke, 1981). 
Together with more accurate bathymetric information the potential of improving 
the analysis of refraction and shoaling in shallow water can be obtained. 


TABLE 8.1: Deviation of INSAR Image Spectra Wavelength and Direc¬ 
tion from Wave Array and Model Prediction Output 



Wave Direction 

Wavelength 


T = 9.1 sec. 

T =s 15.9 sec. 

T as 9.1 sec. 

T = 15.9 sec. 

Correlation 

0.43 

0.88 

0.96 

0.99 

rms error 

6.2° 

L _ _ _ 

6.4° 

3.9 m 

j 

4.7 m 






















m 



Figure 8.5: Predicted refraction versus INSAR 3D wavenumber spectra 
angles for the two observed swells (9.1 sec., 15.9 sec.). 
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INSAR WAVELENGTH (M) 


Figure 8.6: Predicted shoaling versus INSAR 2D wavenumber spectra 
angles for the wave array in situ observed swells (9.1 sec., 15.9 sec.). 
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2. Effects of Airborne SAR System on the Study of Ocean Surface 
Waves 

Since the INSAR technique is based on SAR imagery, some key fun¬ 
damental shortcomings of SAR in imaging ocean surface waves were described in 
Chapter III. D. Through examples from the acquired remote sensing data the effect 
of these shortcomings on INSAR wave imagery is discussed in this section. The 
known SAR coherent motion effects include velocity bunching, acceleration defo¬ 
cus and scene coherence time limitation. The noncoherent scene motion effects, 
which are only important for relatively slow scanning platforms like aircraft, include 
scanning distortion and smearing (look misregistration for multi-look processing), 
a. Scanning Distortion 

Scanning distortion arises due to the relative motion between the 
projected translation of the mean pattern of image reflectivity by the wave phase 
velocity C, and the imaging radar platform velocity V. Therefore, this relative 
velocity effects during the observation time is in general proportional to the ratio 
C/V. 

The image patterns should be corrected if the projected azimuthal 
component of the phase velocity is significant compared to the platform velocity. In 
this case the image contrast can be improved by adjusting the processor focus by 
C*/2. The spectral energies have to be corrected in the azimuthal direction by kC/V. 
Using this method suggested by Vachon et al. (1988), the inherent 180° ambiguity 
(for all imaging radars) in the wave peak directions can be resolved. Artifact spectral 
contours obtained from opposing flight directions will diverge after this correction, 
while the actual contours will converge. Using Equation (3.16) and the computed 
INSAR wavelength and direction from spectral analysis (see Table 7.5), the scanning 
distortion effect is estimated. Longer waves with higher phase velocity are relatively 
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more distorted. Therefore, the longest swell (310 m) computed from station 4 (360° 
image, Table 7.5) was selected to demonstrate the scanning distortion effect. The 
calculated wavelength distortion is of the order of 2% and 5° in direction. For 
shorter waves observed in 360° image where the waves were mainly range oriented, 
the scanning distortion for both.swells is practically negligible. Figure 8.7 illustrates 
INSAR image spectra computed from sub-areas in the wave array and in deeper 
water centered 2000 m west from the wave array zone for the eastern 90° flight 
path. Shoaling and refraction effects are observed for the shorter swell only (peak 
PI) since the long crested swell is less, much less energetic in deeper water and was 
not observed. In shallow water the long crested swell becomes more energetic due 
to shoaling and a relatively weak peak with a wavelength of 212 m and direction of 
74°, 14° norther from the computed (direction 260°) wave array directional spectra, 
was measured. Due to the azimuthal orientation of the swells in the 90° flight path 
the wavelengths are distorted by almost 6-7%, which is relatively a large effect 
compared with the corresponding 360° flight path image in which these swells had 
a range orientation. Correction for scanning distortion will make these swells more 
comparable to the wave array wavelength estimated (192 m, Table 7.5). 

In deeper water a relatively large rotation of the shorter swell to¬ 
wards range direction is observed. This additional rotation (also of the long swell 
in shallow water) is not well understood at the moment. The noncoherent scene 
motion effects are negligible for spacecraft since the ratio C/V is small, 
b. Velocity Bunching 

For azimuthaly oriented waves, the velocity bunching imaging mech¬ 
anism is significant. Though the INSAR image formation is not based on this 
mechanism it is formed from complex SAR images, which might be distorted due to 
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Figure 8.7: INSAR Eastern 90° flight path image spectra computed from: 
Top: Shallow water sub-area close to the Marina wave array; Bottom: 
Deeper sub-area centered 2000 m west from the shallow sub-area. 
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bunching effects. For INSAR imagery it is felt that velocity bunching might con¬ 
taminate the image formation. This imaging mechanism becomes less significant as 
the wave direction become more range oriented, the wavelength increases and/or 
the wave height decreases. Based on the theory described in Chapter III. D and 
INSAR orbital velocity estimates, the velocity bunching effect is evaluated. 

The fractional amplitude of the azimuthal image shift is estimated 
using Equation 3.13. The wave amplitude is computed from the INSAR image 
data, using linear wave theory and simple geometry considerations. The mean 
estimated orbital velocity obtained from 360° INSAR image for the 9.1 sec. swell 
varies between 20-23 cm/sec. This velocity is consistent with the rough estimate 
of Uob using Equation (7.13). The shorter swell was selected for a test case since 
it has a relatively larger orbital velocity and hence experience a stronger bunching 
effect. The imaged area is considered as intermediate water depth, therefore the 
orbital wave particle velocity pattern is elliptic rather than circular. The equation 
of this pattern is given by: 

U 2 W 2 

A* + A^ufikh = 1 ^ 


where A' = AgT'/L , and U, and W are the horizontal and vertical velocity compo¬ 
nents. Transforming Equation (8.7) into polar coordinates after simple manipulation 
the wave amplitude is obtained 


Lr ( sin 2 7 

A= W (p^ + C0S7 


) 


1/2 


( 8 . 8 ) 


with 


Uob 

(sin 2 7 sin 2 (f> + cos 2 <f>) 1 ! 2 


(8.9) 


where 7 is the incidence angle and <j> is the angle between the wave direction and 
the radar platform attitude. For both swells estimated from INSAR image spectra 
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close to the Marina wave array, C^, and C' are computed using Equations ( 8 . 8 ), 
(8.9), (3.13), and (3.14). Table 8.1 shows that both swells are far from being crit¬ 
ical. The shorter swell is more influenced by velocity bunching. The highest ratio 
of C'IC'c,. = 0.3 is obtained for the eastward 90° flight path due to an increased 
azimuthal orientation of the imaged swell. This might be one of the reasons for the 
more smeared INSAR image obtained from the eastward flight path relative to the 
northward 360° flight path image. The results depicted in Table 8.1 mainly show 
that the imaged scattering elements redistribution is one to one resulting in single 
unambiguous brightness patterns in the SAR imagery representing the actual ocean 
waves. 

TABLE 8.2: Actual Velocity Bunching Parameters (C, Co) Estimates Ob¬ 
tained from INSAR Imagery (Northern Flight Path 360°) 


Wave 
Period V 
(sec) 

Flight 

Path 

7 


h 

(m) 

| 

Uqb 

(cm/sec) 

A 

(m) 

O'er 

C' 

9.1 

360° 

40° 

69° 

16 

97 

23 

0.3 

43 

10 

15.9 

360° 

40° 

84° 

16 

192 

15 

0.2 

292 

0.4 

9.1 

90° 

45° 

12 ° 

16 

120 

16 

0.23 

19 

6 


c. Coherent Motion Effect 

The degradation in azimuthal resolution due to finite scene coher¬ 
ence time has a low pass filter effect on the SAR image. It shifts the azimuthal wave 
number to lower values. An example of this effect is inferred by comparing the wave 
array spectra (Figure 7.11, Bottom) and the corresponding INSAR image derived 
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spectra (Figure 8.7, Top). The in situ wave array spectrum clearly shows another 
peak with a wavelength of almost 70 m and direction of 300°. For flight path 90° 
this relatively short swell is basically azimuthally oriented and is not seen in the 
INSAR image spectra. It is concluded that the wave patterns in the INSAR images 
obtained from flight paths 90° and 270° are less pronounced mostly due to the az¬ 
imuthal smear caused by the finite scene coherence time limitation. An evidence of 
this conclusion can also be observed in Figure 7.4. At low incidence angles where 
the degradation in azimuth resolution is relatively small waves are clearly observed 
in the INSAR image. At higher incidence angles where the degradation is larger, 
wave patterns are hardly seen. This evidence illustrates the importance of the finite 
scene coherence time in determining the SAR and INSAR capability to image ocean 
surface waves. 

3. INSAR Derived Spectra Conclusions 

• The experimental results show how interferometric SAR may be used to image 
ocean surface waves. Given the problems associated with imaging waves using 
conventional SAR, this unique technique represents a significant advance. The 
agreement between INSAR image spectra and in situ derived spectra (Figures 
7.11, 7.14) is encouraging. In spite of this potential, this section illustrated 
some of the problems that might arise with an airborne INSAR in imaging 
ocean surface waves. 

• Noncoherent scene motion effects are dependent upon C/V. For a relatively 
slow airborne SAR (and INSAR) platform , noncoherent motion effects might 
lead to wave imaging complications that include scanning smear and distor¬ 
tion, which might result in ambiguous image wave patterns. This effect is 
neglibile for a fast moving platforms like spacecraft. Focusing adjustment and 
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scanning distortion correction methods are not reliable for multimodal wave 
fields, which is usually the actual situation of the ocean surface (Vachon et al., 
1988). 

• The finite scene coherence time degrades the azimuth resolution which is pro¬ 
portional to R/V and was shown to have an important role for ocean wave 
imaging for aircraft and spacecraft. Another coherent motion effect is the ve¬ 
locity bunching. The observed wave field (flight path 360°) contained mainly 
range-oriented waves. This orientation is favorable for radar imagery since it 
minimizes image distortion due to azimuth image shift (velocity bunching). 
The imaged waves in the Marina experiment did not experience significant 
coherent azimuth image shift (Table 8.2). 

• Several effects degrade the INSAR image wave patterns (see Figure 7.4). The 
dominant effects are believed to be: the degradation in azimuth resolution 
due to finite scene coherence time, looking crosswinds rather than upwind 
and speckle and random noise. The lesser effects for the Marina experiment 
conditions were scanning distortion and probably velocity bunching. 

• In shoaling waters the area of the image being subjected to spectral analy¬ 
sis should be reduced to allow the waves to be represented locally without 
refraction and shoaling effects. 

E. SOURCE OF SWELLS 

The Global Spectral Ocean Wave Model (GSOWM) outputs for August 31 
and September 8, 1989 dearly illustrate the possible sources and directions of dom¬ 
inant wave propagation (see Figures 6.8-6.11). A GSOWM 2D spectrum taken at 
a grid point about 60 miles southwest off Monterey Bay (Figure 8.8) indicates a 
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9.5 second northwesterly wave system together with a narrow-band southerly swell 
(16 sec.). The propagation of the long crested southerly swell into Monterey Bay 
is consistent with the increasing importance of refractive features with decreasing 
wave frequencies in the southern edge of Monterey Bay. 



Figure 8.8: GSOWM directional frequency spectrum obtained from a grid 
point at 35° N Lat and 122.5° W Long on September 8,1989. A southerly 
16 sec. swell and a 9.5 sec. northwest wave system are observed. 


The characteristics of southerly swells in the northern hemisphere are not 
well understood. Munk et al., in a series of papers, established the source as high 
latitude storms in the South Pacific and Indian Oceans. These swells are usually less 
energetic (smaller amplitude after traveling 10,000 km) and longer in period than 
swells originating in the North Pacific (Guza, 1986 and Semdin, 1986). Based on the 
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global models of pressure fields, winds and spectral analysis, it is assumed that the 
south-southwest long crested swell observed in the Marina experiment, September 
8 , 1989, was generated by a high latitude southern hemisphere storm and that the 
west-northwest shorter period swell was generated by a northwest Pacific storm. 






IX. CONCLUSIONS 


The present study reveals the considerable ability of this novel technique 
to provide reliable quantative spatial information about dynamic processes at the 
ocean surface. Realization of the potential of Interferometric SAR with spacecraft 
might make it possible to monitor global ocean wave spectra and surface currents. 

The mean estimate of the L-band scene coherence time 0 (100 msec) obtained 
from INSAR and SAR imagery during the Marina experiment, on September 8,1989, 
is consistent with the sparse known observational estimates obtained under similar 
moderate sea state and ambient conditions. This parameter is shown to be the 
dominant limiting factor in determining the SAR and INSAR azimuth resolution 
and hence the azimuthal cutoff wavenumber for imaging ocean surface dominant 
waves. Waves with azimuthal component shorter than 0 (100 m) could not be 
observed. The azimuth image shift (velocity bunching) caused by the coherently 
sensed radial component of the orbital velocity of the long waves might disturb the 
INSAR phase measurement. It is shown that the noncoherent scanning distortion 
effect is important if the projected phase velocity of the imaged wave is significant 
relative to the platform velocity. This effect is negligible for fast scanning platforms 
like spacecraft. It is concluded that range oriented waves are favorable for INSAR 
(and SAR) where the undesired motion effects like the temporal decorrelation of the 
ocean surface and velocity bunching do not distort the image. Also, range waves 
project a larger radial component of orbital velocity (which is sensed by INSAR), 
compared with azimuthal waves. Additional theoretical and experimental studies 
are required to quantify and to provide better understanding of these issues in terms 
of INSAR performance to image ocean surface wave fields. 
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The results of this experiment in estimating ocean surface currents, together 
with earlier work by Goldstein et al. (1989), and for the first time ocean wavenumber 
spectra suggest some additional conclusions: 

• The method is shown to be most reliable when the radar observation is along 
the direction of the strongest current component. 

• The interferometric technique appears to be more accurate when the radar is 
looking upwind (and probably downwind). In these cases, the phase velocity 
bias from the Bragg waves responsible for the radar signal is less ambiguous, 
and the returned amplitude signal is strongest. 

• Shoaling and refraction effects are clearly observed in the INSAR phase images. 
In shoaling waters, the area of the image being interpreted via spectral analysis 
should be reduced to insure homogeneous waves. 

• The interferograms contain numerous additional ocean surface features which 
were not pursued in the present study, such as convergence, divergence, and 
vorticity fields and internal and other long waves. 
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APPENDIX A 

INTERFEROGRAM PROCESSING 


The recorded SAR data from each antenna consists of a set of four files. Each 
file consists of 1024 azimuth by 750 range pixels. P m y and Q m ij are the corre¬ 
sponding complex compressed pixel values in the imaging map as processed from 
the recordred data received from each of the two antennas. The subscripts m,i,j 
are the file number, row, and column of a pixel, respectively, 7y is the complex 
image value of a pixel in the interferogram. Ay, $y is the magnitude and phase of 


Processed Imaging Files from each Antenna 




/y is derived by taking the complex conjugate product of the sum of the four 
corresponding original complex image pixels obtained from each antenna. 
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m=l 

(A.1) 

Aij exp[t$,j] 

(A.2) 


The interferogram file as received from JPL is processed to a 2048 azimuth 
by 750 range pixels. The first 1024 azimuth pixels are magnitudes of the complex 
reflectivity (Ay) and the other 1024 pixels are phases which consists of information 
of surface velocity. 

The pixel averaging in Equation A.l is in expanse of azimuth resolution. The 
pixel size is the resulting interferogram is 12.1 m is azimuth and 6.6 m in the viewing 
direction. 
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APPENDIX B 

ESTIMATION OF RANDOM ERROR IN 
LAGRANGIAN DRIFTER VELOCITY 


A. THEORETICAL BACKGROUND 

Drifter velocity is a simple function of distance traveled divided by time, 
where distance is found by differencing positions. Thus, random errors in the posi¬ 
tions do propogate into the velocity. 

The velocity, V, is a function of the starting and ending position, Xi and X 2 , 
and the time difference AT = T 2 — Tj 


V = F(X 1 ,X 2 ,AT) (BA) 

where V, Xi and X 2 are all vector quantities with east and north components given 
by: 

V = (u,v) X J = (x 1 ,y 1 ) X 2 = (x 2 ,y 2 ) (B.2) 


Specifically, 

u= (x 2 - x x )/AT v = (y 2 -y x )/AT (B. 3) 

The error propagation equation can be set up with £(V) being the 2x2 
variance-covariance matrix of the velocity components 

£(V) = G £(X„Xj) G t = [ °* I (BA) 


where G is the 2x4 Jacobian of partial derivatives and is defined as 


G, 


dFi(u,v) 


* = 1,2 




55 j (X i ,X 2 ) j = 1,2,3,4 


or simply in our case 


G = 


10-1 0 
0 1 0-1 


1 

AT 


(B. 5) 


(B.6) 
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and E x is the 4x4 combination of variance-covariance matrix for the variable 
positions and is taken to be 


S(X„X 2 ) 




0 

0 ' 


< 

0 

0 

0 

0 

a lx 

O'*! VI 

0 

0 

<TxiVi 

< . 


(B. 7) 


B. OBSERVATION ERROR ESTIMATES 

Using Equations (B.l) through (B.7), the velocity error estimates are com¬ 
puted: 

1 . Group I: Between observation 2 and 1 


T\ = 11,49,30 (hours, minutes, seconds) PDT 
T 2 = 13,15,44 PDT 


A T = 
S(X.) = 

E(X,) = 

£(V) = 

<7u = 


86.2 min 

20.07 2.69 
2.69 8.24 

22.76 2.69 
2.69 S.03 

1.6-lO" 2 2 • 10 -3 
2 • 10 -3 6 • 10 -3 

0.13 cm/sec 


m 


m 


cm sec 


<j v — 8 • 10 -2 cm/sec 


2 . Group I: Between observations 3 and 2 

T z = 14,16,12 PDT 
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T 2 = 13,15,44 PDT 


= 60.5 min 


S(Xi) = 


S(X 2 ) = 


E(V) = 


25.59 2.7 
2.7 7.88 

22.76 2.69 
2.69 8.03 


3.7 • 10~ 3 4 • 10- 3 

4 • 10“ 3 1.2-IQ" 3 


= 0.19 cm/sec 
a v = 0.11 cm/sec 

3. Group II: Between observations 2 and 1 

T\ = 12,05,28 PDT 


T 2 = 13,36,08 PDT 


= 90.7 min 


S(X x ) = 


E(Xt) 


E(V) 


27.2 2.13 

2.13 7.67 


32.26 1.49 

m‘ 

1.49 7.44 

2 • 10" 2 10~ 3 
10" 3 5•IQ -3 


2_„_—2 


cm sec 


= 0.14 cm/sec 
= 7.10 -2 cm/sec 


4. Group II: Between observations 3 and 2 


2 -2 
cm sec 


T 3 = 14,42,39 PDT 


180 


T 2 = 13,36,08 PDT 


AT = 6.65 min 


E(X.) = 

33.78 0.24 

0.24 7.33 

m 2 

E(Xj) = 

‘32.26 1.49 ' 

1.49 7.44 

m 2 

S(V) = 

' 4.1 • 10- 2 10- 3 

10~ 3 9 • 10 -3 

2 -2 
cm sec 

cr u = 0.2 cm/sec 


a v = 0.1 cm/sec 
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APPENDIX C 


IMAGE COORDINATES GROUND 
TRUTHING USING GPS DATA 

A. IMAGE ADJUSTMENT PREPARATION 

A modified Transverse Mercator projection was used to convert the latitude 
and longitude ground coordinates to a plane coordinate grid of eastings and nor¬ 
things, ( E,N ) coordinate system (Floyd, 1985). The projection was centered in 
the image area in order to minimize grid distortions. GRS 80 ellipsoid values were 
used which are essentially the same as for the WGS 84 geometric ellipsoid. The 
projection parameters are: 


Ellipsoid Semi-Major Axis, 
an Inverse Flattening, 1/f 

6378137.0 meters 
298.257222101 

Central meridian (longitude) 
False latitude 

False easting 

Scale factor, central meridian 

W 121° 47.0’ 

N 36° 40.0’ 

3000 meters 
1.00000 


The range pixel number ( y ) of the image control points needed a small cor¬ 
rection for elevation distortion (Ulaby, 1982). The corrected range pixels, y c , were 
computed as follows: 

y lc = 1 IK (^ v (y- l) + r t> )2 + 2/fs-s2 - r a ) + 1 (C.l) 

where 
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TABLE C.l: Ground and Image Coordinate Data 








Type 

Latitude (d m s) 1 
Longitude (d m s) 
Elevation (m) 

Plane 

Coordinates 

(m) 

Easting 

Northing 

Image 

Coordinates (pixels) 

X 

y 

A90 A270 A360 

Bridge 

N 36 43 58.34 

3204.3 


313 

197 

abuttment 

W 121 46 51.77 
6.1 

7347.1 


139 

69 

Road 

N 36 43 15.46 

2630.4 

146 

295 

245 

intersection 

W 121 47 14.89 
15.2 

6025.0 

484 

222 

104 

Fence 

N 36 42 47.87 

4169.5 

82 

360 


corner 

W 121 46 12.88 
30.5 

5174.9 

434 

280 


Fence 

N 36 42 37.74 

3645.5 

104 

339 

296 

corner 

Fence 

corner I 

W 121 46 33.99 
33.5 

N 36 42 12.89 
W 121 46 16.03 
42.7 

4862.4 

4091.3 

4096.4 

415 

86 

364 

305 

24 

Road 

N 36 41 15.95 

1568.8 

190 

250 

397 

intersection 

W 121 47 57.65 
3.0 

2341.3 

261 

493 
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IMAGE 

PARAMETERS 

Azimuth Scale, k x \ 

(meters/pixel) 

24.20 

24.20 

24.20 

(a priori values) 

Range Scale, k v (meters/pixel) 

16.14 

14.54 

12.34 


Flight azimuth, 

a (degrees) 

90° 

- 

270° 

000° 

IMAGE 

CONSTANTS 

Flying height, H (meters) 

8285 

8273 

8178 


Ground range to 1st pixel, r 0 
from nadir (meters) 

2601 

2918 

4015 


1 d m s = deg min sec 
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y = range pixel number (uncorrected) 
k v — a priori range scale (meters/pixel) 
r 0 = ground range from nadir to y — 1 pixel (meters) 
h = flight altitude (meters) 
z = elevation of ground point (meters) 

The magnitude of these corrections was small, generally less than 2 pixels. However, 
for a few points at elevations of 30 meters or more near the flight line, the corrections 
ranged from 3 to 5 pixels. The correction always increased the pixel number, moving 
the point further from the flight line. Only the control points needed this correction, 
the imaged water surface assumed to .be at zero elevation were not subjected to 
elevation distortion. 

B. IMAGE ADJUSTMENT COMPUTATION 

The adjustment can be viewed as solving for the best mathematical mapping 
between the image and ground coordinate systems, between (x,y) and (E, N), sub¬ 
ject to assumptions about the underlying form of the mathematical transformation. 

The adjustment used a least-squares technique in which the observed data are 
weighted inversely proportional to their variances, as described in most surveying 
or photogrammetry texts. This weighting strategy results in a minimum variance 
solution, which is generally considered the best solution when using data of unequal 
precision. The specific form applied here is called the “combined observation and 
condition equation” method (Uotila, 1986). 

The assumed transformation was a simple scaling, rotation, and translation 
between the (x,y) and ( E,N) coordinates systems, with the addition of separate 
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scales in the x and y directions (i.e., the pixels were assumed to be rectangular, not 
square). The assumed form of the transformation equations was: 


Ei - E 0 + Xik x sin a + yik v cos a — 0 (C.2) 

Nj - N 0 + Xik x cos a - y{k v since = 0 (C.3) 


where 


Ei,N{ = ground coordinates of ith control point, 
x,-, yi = image coordinates of ith control point, 
k y — scales (meters/pixel) in x and y directions, 
a = azimuth of the flight line and rotation between 
coordinate systems(a = 0 points to north), 

E 0 , N 0 =s ground coordinates of image origin, (x,y) = (0,0). 


Note that the x-axis points are opposite to the flight line direction because of the 
way the image processing hardware loaded the image data. This accounted for the 
unusual appearance of the signs in these equations (C.2, C.3). Also note that a is 
measured clockwise from north, which causes the role of the sine and cosine factors 
to be reversed in comparison with some other presentations. 

In the adjustment terminology, the quantities k x ,k y ,a,E 0 , and N 0 are pa¬ 
rameters to be solved for. Hence there were 5 parameters in each image. For each 
control point, one pair of condition equations (C.2 and C.3) were formed. As long 
as the number of condition equations exceeded the number of parameters, the least 
squares solution was possible. 

Three addition condition equations were used, reflecting the a priori estimates 
of the parameters k x , k y , and a: 


k x k X o — 0 


(C.4) 
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(C.5) 

(C.6) 


ky kyo ~~ 0 

a - a 0 = 0 

where k xo , k yo , and a 0 were the a priori values of the respective parameters. 

Five control points were used on images A90 and A270, and four control 
points on image A360. Hence, 13 condition equations were used for A090 and R270, 
and 11 for A360; more than enough for the least squares adjustment. Weighting 
required estimation of the standard errors (precision) on each observation. The 
estimated values were as follows: 


<te ~ — 30 meters 

o x — c y = 3 pixels 
cr a0 = 1 degree 
< 7 * JO = 0.1 meters/pixel 
a ky 0 = 0.2 meters/pixel 

The adjust proceeded as follows. Define the set of condition equations as F(L a , X a ) = 
0, where the vector L a represents the adjusted values of the observables AT*, £,•, 
Vi> k xo , k yo) a 0 ) and the vector X a represents the adjusted values of the parameters 
(k x , k y , a, E 0 , N 0 ). The set of condition equations F(L a , X a ) has the form of 
Equations (C.2) through (C.6). 

With, 

N = number of control points on the image 
R = number of condition equations 
U = number of parameters = 5 


then, 


R = 2N + 3 


(C.7) 
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which is the number of condition equations in the set F(L a ,X a ). The length of L 
is 4 N -f 3. The degrees of freedom in the adjustment is R - [/, 

dof = 2N + Z — 5 — 2N-2 = 2(N - 1) (C.8) 

Because the set F(L 01 X 0 ) ^ 0 when L 0 is the raw observations and X 0 is the 
approximate values of the parameters, define a non-zero vector W, length R, 

W = F(L 0 ,X 0 ). 

The adjusted parameters are found by iterating the following matrix equation until 
the X a vector converges, 

X a - X 0 + dX (C.9) 

dX = -{A'M-'AY'A'M-'W (C.10) 

where matrix A = dF/dX (R rows by U columns), matrix M = BCB 1 (R by i?), 
matrix B = dF/dL (R by L), and matrix C = covariance matrix of observables (L 
by L). 

In this case, C was a simple diagonal matrix of estimated variances of each 
observed quantity, and provides the relative weighting of the observations in the 
adjustment. The traditional weight matrix is called P, where P = C~ l . 

For analysis of the quality of the adjustment, residuals on the observed quan¬ 
tities were computed, as a vector 

V = L a -L 0 = -CB t M- 1 (AdX a + W) (C.ll) 

and the “ variance of unit weight ” was computed as the scaler 


- (V'PV)ldof. 


(C.12) 




The adjustment also computed precisions for the derived values of the parameters, 
expressed as a covariance matrix of the parameters, C x , 

C x = (A* ATM) -1 (C.13) 

Each of these adjustment results, Equations (C.9) through (C.13), were com¬ 
puted for each image. The results are summarized in Table C.2. 

C. POST-ADJUSTMENT IMAGE SHIFT 

Normally, an adjustment process like that described would be sufficient to 
accurately align an image with ground coordinates. In this case, however, the ad¬ 
justed results might have been degraded by the lack of available control points. In 
typical image adjustment problems the control points are selected to span the en¬ 
tire image area. Because most of each image was over water, the control points 
were concentrated in a relatively small land portion. Hence, the danger the results 
might not be reliable in the water area. Essentially, the danger was one of using 
extrapolation rather than interpolation. 

Examination showed that the adjusted scale and azimuth values agreed rea¬ 
sonably well with their a priori values — but not exactly (most within 1 standard 
deviation). But because these three parameters determined the overall size and ori¬ 
entation of the image, it was decided to accept the a priori values rather than the 
adjusted values. The a priori values were derived from flight navigation data and 
known characteristics of the SAR instrument. Consequently, it was felt they better 
represented the true values for the enthe image area. 

This left the adjustment to determine the location of the images. The im¬ 
age origin coordinates ( E 0 ,N 0 ) solved for in the adjustment reflected the adjusted 
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TABLE C.2: Image Adjustment Results 



Aljj ky 

a 

E 0 

N 0 

Image 

Parameter (m/pix) (m/pix) 

(^g) 

(m) 

(m) 

A90 

Initial 24.20 16.14 

90.0 




Adjusted 24.20 16.20 

89.9 

6173 

-1864 


a (0.04) (0.07) 

(0.3) 

(35) 

(33) 


Post adjustment origin shift 


-8 

5 


Final origin coordinates 


6165 

-1859 


Variance of unit weight 0.14 





Latitude limits N 36-39-00 to N 36-43-28 




Longitude Limits W 121-44-54 to W 121-53-12 


A360 

Initial 24.20 12.34 

000.0 




Adjusted 24.23 12.51 

001.3 

4129 

11988 


a (0.12) (0.22) 

(0.8) 

(91) 

(60) 


Post adjustment origin shift 


-384 

68 


Final origin coordinates 


3745 

12056 


Variance of unit weight 1.76 





Latitude limits N 36-39-49 to N 36-46-30 




Longitude limits W 121-46-30 to W 121-50-44 
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TABLE C.2: Continued 


Image 

Parameter 

k x 

(m/pix) 

ky 

(m/pix) 

a 

(deg) 

E 0 

(m) 

N 0 

(m) 

A270 

Initial 

24.20 

14.54 

270.0 

_ 

_ 


Adjusted 

24.19 

14.34 

271.0 

-4416 

9437 


a 

(0.09) 

(0.14) 

(0.7) 

(68) 

(105) 


Post Adjustment origin shift —69 —58 

Final Origin coordinates —4485 9379 

Variance of unit weight 1.39 


Latitude limits N 36-41-03 to N 36-41-03 
Longitude limits W 121-43-42 to W 121-52-01 

rather than a priori image size and orientation. The further the origin point from the 
actual control points used, the less precise its ground coordinates were determined. 
To compensate for this effect, ground coordinates of the four image corners were 
computed in two ways, using (1) the adjusted and then (2) the a priori parameter 
values. The origin coordinates were then shifted to align the centers of the two 
versions. In other words, any systematic errors were forced to the edges of the 
image and minimized near the center. 

The amount of this post-adjustment origin shift was barely significant on im¬ 
ages A90 and A270, in comparison with the standard error of the adjusted coordinate 
values (see Table C.2). 

On image A360, however, this shift was large, about 180 meters westward, 
which was twice the estimated standard error. No obvious explanation for this 
discrepancy was found, and subsequent comparisons of the predicted wave gage po¬ 
sition with the distance from the shoreline in the image indicated that an additional 
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200 meter westward shift was needed; for a total shift of about 380 meters west of 
the adjustment derived position. Whatever the error, it did not appear to affect the 
north-south direction. 
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